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Chapter 1 



The Turbulent Refractive Index: 
Dynamics and Stochastic 
Properties 

The study of phenomena occurring in a turbulent fluid has been successively improved 
during the last 40 years. Specifically, the concentration of a substance advected by the 
turbulence has received most of the attention, for it covers a wide range of natural and 
engineering settings: heat transport, dye diffusion, microscopic organism movements, 
etc.. These substances are described by scalar fields with a negligible back-effect on 
the flow; thus, they are called passive scalar fields. 

The turbulent refractive index also belongs to this class; this is not a novelty 
(Tatarski, 1961). The temperature is a passive scalar held whenever it produces buoy- 
ancy forces smaller than the inertial stresses driving the flow, and a direct calculation 
shows that its fluctuations are proportional to those of the index. 

Our interest in lightwave propagation through turbulent media must start here 
then. That is, we have to comprehend the media before attempt a description of 
the propagation itself. In the forthcoming sections we will study the dynamics and 
stochastic properties of passive scalars, and eventually propose models for the refractive 
index. 
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1.1 Turbulence 

1.1.1 The turbulent flow: Kolmogorov hypotheses 

Above wc have, without more precisions, referred to the turbulent media. From now 
on, we mean incompressible fluids in turbulent state; moreover, all our discussions will 
be targeting the atmosphere or experiments that resemble it. Of course, this section is 
intent to explain what Hurbulenf is. 

Let us start from the beginning; as it is well known fluids are governed by the 
Navier-Stokes equation: 

d 1 

— v+(v- V)v-i/Av= -(F- Vp), (1.1) 
at p 

v(r, t) -.M? ^ M_|_ — > M? is the velocity fleld, while v is the viscosity of the fluid (with 
dimensions [v] = L'^/T), p the density, p the pressure and F the external force. It is 
worth noting that this equation is scale invariant. So it can be turned into the following 
adimensional equation, 

|,+ (,.V)v-(^)Av = i(F-Vp), (1.2) 

with I and r the characteristic length and time of the system. The constant multiplying 
the first term at the right-hand side of the latter equation introduces the Reynolds 
number, 

7^e(/) = — , (1.3) 

V 

vi is the velocity change on the scale length /. The Reynolds number is a scale dependent 
quantity, and its magnitude measures the flow regime: it compares the non-linear 
advection term (v • V)v against the dissipation —u Av. While low Reynolds numbers, 
TZe{^) Ij correspond to regular and laminar flows and intermediate numbers, 1 < 
TZe{l) ^ 10^, exhibit complex patterns, higher Reynolds numbers, TZeiO ^ lO"'; drive 
the flow to an apparent spatial disorder: parcels of fluids follow chaotic trajectories. 
In particular, when the Reynolds number tends to inflnity the flow exhibits a fully 
developed turbulence. The non-linear advection is preponderant because the dissipative 
term goes to zero. 
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The equation (1.1) induces the energy balance equation (per unit mass): 




The balance is given here between the first term on the rightmost-hand side of this 
equation, which represents the energy injected per unit time into the system, and the 
energy dissipated by the viscous forces, that is, the second term. 

It was Kolmogorov (1941) who first reahzed that from dimensional and rehable 
heuristic arguments the energy transfer could be explained. His success was to notice 
that the results of this analysis become universal laws in the statistical sense. The 
turbulent velocity field should be thought a stochastic variable in the ensemble's sense 
of the statistical mechanics. It is independent on how the turbulence began: it does 
not matter the way the energy is injected. That is, the statistics of the chosen force 
has no effect over the statistics of the turbulence. 

Moreover, we will also assume that a fully developed turbulence is spatially isotropic, 
homogeneous, and stationary: for any linear transformation and translation the system 
looks the same. 

In this section we will treat the turbulence development under the direct energy 
injection. That is, the energy is injected by the largest disturbances of size L — the 
integral scale — , corresponding to the size of the bath, and then it is transferred towards 
the smallest scales. Finally a minimum scale Iq — the inner length — is reached, there 
the energy is dumped by the viscosity into heat (the magnitude of the inner length 
oscilates between lO^^m and 10~^m). 

The range of scales I where the energy transfer happens without loss, the flux of 
energy from scale to scale is constant, is called inertial range 

This process can be thought as a cascade of energy that propagates through the scales 
via a succession of disturbances {eddies which are portions of fluid with size / and 
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velocity vi): big eddies break up smaller ones. These eddies are arranged in a hierarchy 
according to its size, from the bigger to the smallest, as follows: 

Z„~L7r", n = 0,l,--- (1.5) 

with TT < 1 the contraction ratio of the eddy size from one generation to the other. 

Now, we can use this scheme to estimate some of the quantities involved in the 
generation of the turbulence. Thus, let Vi^ be the volume occupied by the eddies of 
the n-th generation; their energy density is ui^ ~ ''^i„/'^- straightforward then that 
the accumulated total energy by the eddies of size I ~ Z„ is. 

El - vfVi. (1.6) 

Knowing that the characteristic life-span of the disturbances is r; ~ In/vi- We obtain 
the following estimation for the energy transfer rate. 

El vfVi 

£ ~ ~ ; . 

n I 

If we now consider that the volume occupied by the eddies is independent of the scale, 
i.e., Vi ~ const.. Then the energy flux per unit volume e is also constant and, 

el^vf. (1.7) 

This scaling law is the fundamental result, as we shall see, of the whole chapter for it 
will be underneath every property we are about to show. 

For instance, let us try some examples: we have defined before the inner length as 
the scale where the dissipative term becomes noticeable with respect to the convective 
one, that is, 

II (V • V)v|| ~ VI, ~ e'/' I-'/' ~ II - . Av|| e^'- lo"\ 

and thus the inner scale is roughly 



(1.8) 
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It vanishes as — > 0, and this results to be an ultraviolet cut-off. Below this cut-off 
the advection can be neglected and the velocity turns more regular. 

The local Reynolds number at a scale I can be calculated from (1.3); it is 7le{l) ~ 
{l/hY^^- Moreover, the system's Reynolds number may be taken as TZe :— TZe{L). So, 
the condition I ':$> Iq is in agreement with the conditions for turbulence development: 
the inertial range grows as the system's Reynolds number do so. 

Prom equation (1.7) we can also check the occurrence of equilibrium between the 
injected and dissipated energy. Using the isotropy and homogeneity properties of the 
velocity field: 




Let us start looking at the stochastic properties of the velocity field. That is, we 
want to compute the n-point correlations of the turbulent velocity. The usual procedure 
to overtake is as follows: first, we separate the (stochastic) fiuctuations u from the 
mean (averaged) fiow (v), so it can be written v(r,t) = (v)(r) + u(r,t); second, we 
derive from the Navier-Stokes equations the corresponding equations for the n-point 
correlation. Here the real problem arises: these equations are non-linear, for it is a 
closure problem. Their solutions are found only in approximation. The Kolmogorov's 
method is so successful because it allows us to override this second step. 

Assume the mean fiow is zero and the random velocity field has the properties we 
have discussed at the beginning: homogeneity and isotropy. The existence of scaling 
laws for the n-correlation functions means the existence of exponents rjn such that 

31imlim r''"(u(/ ri, t) • • ■u(/r„,t)). 

Because the energy transfer per unit volume is constant all over the inertial range 
and there is independence from the source of turbulence, the n-correlators of the 
stochastic velocity field will just depend from the scale. If we look back at equation 
(1.7), we have: 

(||u(r + r')-u(r')ir) = C„(£||r||)"/3, (1.10) 

and the constants C„ are universal. 

This scaling method, although effective determining some fundamental behavior 
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of the turbulent velocity field, is scarce explaining the way the transference of energy 
occurs between the different scales. This the task to be tackled in the next section. 

1.1.2 The energy cascade in isotropic turbulence 

As before, we are dealing with an homogeneous and isotropic turbulence. The energy 
is then constant throughout space. Thus, when we consider the transport of turbulent 
energy, this will be in wavenumber rather than in the coordinate space. So, we can 
foresee a transfer from one range of eddy sizes to another: the cascade phenomenon. 
We used McComb's (1991) book as the main guide for this section. 

In order to have isotropy and homogeneity we will make the boundary of our system 
go to infinity. We will deal here and thereafter with incompressible fluids, so going back 
to equation (1.1) we set V • F = 0. Additionally, we can obtain another property from 
the Navier-Stokes equation, applying the divergence to both sides of it yields: 

^(^-V) = - E ai-git = - E = -(V « V) . (u « u). 

j,k j,k 

where (8) is the tensor product. That is, each vector is understood as a column matrix, 
r e R X R^, and the inner product acts column by column hke above. This is a Poisson 
equation, and it can be solved calculating the Green's function: 

AG(r,r')=5(r-r'), (1.11) 

with condition n • VG — > 0, as the boundary goes to infinity. The pressure can be 
written, 

p-V(r, i) = -(V ® V) • / dV ^(r, r') (u(r', t) ® u(r', t)) . (1.12) 

There, the superficial terms are zero according to the conditions imposed to the tur- 
bulence. 

The formation of a stationary isotropic turbulence requires the external force F 
to counter-effect the action of the viscous force, but for the present discussion we 
momentarily set it equal to zero. 
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Using the latter equation it takes some effort turning equation (1.1) into: 

^-z/Au = -(l®V)-|(u®u)-V (1®V)- / (iVG(r,r') (u(r') ® u(r')) 
= - (1 (g) V) • D(l (g) V)(u (g) u), 



here (l)jfc = ^jk- The right-hand side of this equation can be changed into a symmetric 
form with the aid of the operator, 

M(V) = -I [(1 ® V) • D(l ® V) + (V ® 1) • D(V ® 1)] , (1.13) 
so we finally find: 

dw 

— -i/Au = M(V)(u®u). (1.14) 
ot 

This equation concentrates all the non-linear effects producing the advection on its left 
side, while the smoothing diffusive term is on the right-hand side. Also, for all the 
practical problems the non-linear term here is no more complex than the original one. 

Although possible, we would rather not build differential equations for the moments 
of the velocity field from equation (1.14); instead, it will be enough for us to recover 
an energy balance equation. Thus, we introduce the Fourier transform of the random 
velocity, 

u(r, t)= \ u(k, t) exp(ik • r). (1.15) 

The continuity equation for incompressible fluids changes in the wavenumber space to 

k-u = 0, (1.16) 

that is, the wavenumber vector is perpendicular to the velocity field. 

We can transform the Navier-Stokes equation into a wavespace equation, as usual, 
using the Fourier Analysis: 

(I- + = ^(k) • / d^k! u(k') ® u(k - k'), (1.17) 



where 

MCkl = 

2i 



M(k) = — [k D(k) + D(k) (8) kl (1.18) 
2z - -I 
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and 



D(k) = I - 



k® k 



:i.i9) 



here I G M'^ x and we have dropped the time dependence on u to simphfy things. 
These last two operators are, of course, Fourier transforms of its counterparts in (1.13). 
It is straightforward from the equation (1.15) that: k ■ D(k) =0: it is a projection. 

Moreover, the moments of u inherit some properties from the turbulent system 
initial setup. In fact, from 



^(^) = 7^T~^ I "(^) exp(-ik • r), 



it follows that the two-point correlation in the k-space is related to the corresponding 
in the r-space by 

(u(k) ® u(k')) 

= — ^ / rfV(iV(u(r)«)u(r'-r))exp[-i(k + k') •r]exp(ik'-r'). (1.20) 

The space correlation is invariant under translations due to the homogeneity of the 
turbulence; so, we have 



(u(r') ® u(r' - r)) = (u(0) ® u(r)). 



:i.2ii 



Equation (1.20) becomes. 



(u(k) ® u(k')) = / rfV rfV (u(0) ® u(r)) X 



X exp[-i(k + k') • r] exp(ik' • r') 

1 



1 



6{k + k') 

+ L(2.)3 

5(k + k') Q(k), 



d'^r (u(0) (g) u(r)) exp(ik' • r 



M3 



(i^r Q(r) exp(ik' ■ r' 



(1.22) 



Q(r) is the isotropic correlation. Hence, the 2-point spectral correlation has a non- 
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vanishing contribution only when k + k' = 0. 

Also, we can prove, with the same arguments, that higher order correlations have 
the same property. That is, 

(u(ki) (g) u(k2) (g) • • • (8) u(k„)) = unless ki + ka H h k„ = 0. (1.23) 

But, it is the isotropy which provides us with what can change all these tensor forms 
for the moments into 1-dimensional expressions. As we said, we are concerned with the 
energy transfer so just the Fourier transform of the second moment will be considered. 
A 2-tensor invariant under rotations and translations can only be expressed as follows 
(Batchelor, 1971), 

Q(k) = B{k)I + A{k)l!i®k 

where the functions A and B are indeterminated but continous. If we multiply this 
equation by k-, and make use of (1-16) then 

k • Q(k) = = B{k)k + A{k)k'^k = [B{k) + A{k)k^]k 

for all k. So, 

q{k) = B{k) = -eA{k), 

and finally it yields 

Q(k) = q{k)l - ^(k ® k) = B{k)q{k). (1.24) 

Now we can make some considerations about q{k). Because trD(k) = 2 from its 
definition. It is 

trQ(k) = tr[D(k)g(A;)] = 2q{k). 

This trace can also be linked to the energy E per unit mass of fluid. The isotropic 
correlation is naturally related to the density of energy, and it gives the following: 

2E - 3{u^) = trQ(r)|^^o = tr / d^kQ{k) = tr / k'^dkq{k) / dVt D(k), 

Jrs Jo j 

here dVL is the solid angle. We have used definition (1.24), and its Fourier relation with 
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the isotropic correlation. The angle integration can easily be carried out, so 

A rco poo poo 

E^—{trI) / edkq{k)^ / ink^dkq{k)^ / dkE{k). (1.25) 
3 Jo Jo Jo 

We have thus defined E{k), the wavenumber spectrum, as the contribution to the total 
energy from harmonic components with wavevectors lying between k and k + dk. The 
quantity q{k) is the density of contributions in wavenumber space to the total energy; 
we will call it spectral density. 

It is now time to calculate the dynamics of the spectral correlation. We will consider 
single-time moments. Henceforth, we (8)-multiply equation (1.17) by u(— k, t) to build 
the matrix. 



_ u(-k, t) ) + z/r (u(k, t) ® u(-k, t)) = 



= M(k) • / d^k' (u(k') ® u(k - k') u(-k, t)). 
We find a similar equation when the change k — > — k is made, that is. 



u(k, t) (8) — ) + ^k'{u{li, t) ® u(-k, t)) = 



= [ ci3yt'(u(k,i)(8)u(k')<8)u(-k-k')) •M(-k). 
Summing both equations, and using the property 

u(-k,t)) + (u(k,^) ^ ^^^) = |(u(k,t) u(-k,t)), 

we finally have: 

(|- + 2i/A;M Q(k,t) =M(k) • / d^k'Q3{k',k-k',t) 

+ / d^fc' Q3(k', -k - k') • M(-k), (1.26) 

Jr3 

here we have defined the 3-point spectral correlation Q3(k',k — k', — k) = (u(k') ® 
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u(k - k') (8) u(-k, t)) e X ]R3 X W. Taking the trace operator and multiplying by 
27rA;^ on both sides of (1.26) we arrive to 



^ + 2uk^^ E{k,t) = T{k,t), 



[1.27) 



where the non-linear term on the right is given by, 
T(k,t) =27rA;2tr jM(k) • / d^k' x 



X 



Q3(k', k - k', -k, t) - Q3(k', -k - k', k, t) 



(1.28) 



This term causes the advection of the spectral energy density: it redistributes the 
energy in the wavenumber space. Henceforth, it should satisfy 



Jo 



dkT{k,t) =0. 



(1.29) 



To prove it let us first notice that from (1.16) and definition (1.19) we have, 

D(k) •u(k) = u(k). 

This property, together with equation (1.16), induces 

M(k) ■ Q3(k', 1, -k) = D(k) ■ Q3(k', 1, -k) = tr i,3Q3(k', 1, -k). 



:i.3o) 



:i.3i) 



where tr 1^3 is the trace computed from the first and third arguments of the 3-tensor 
leaving the second free. 

Now, the integral (1.29) can be put in the form: 



/ dk2T(k,t)^ I Ank^dk I d^k'(2i)-^x 
Jo J Jw^ 

tr i,3Q3(k', 1, -k) - tr i,3Q3(k', 1 - 2k, k) 



xk 



:i.32) 



It is also true 1 ■ tr 1 3Q(1) = 0, for it is 1 ■ u(l) = 0. We can replace k in the first term 
on the right-hand side above by k — 1 = k', because of condition (1.23). Therefore, we 
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find 

poo p p 

/ dkT{k,t)^ / d^k / d^k'{2i)-^x 

Jo JrS Jr3 

X k' • tr i,3Q3(k', 1, -k) - k • tr i,3Q3(k', 1 - 2k, k)] . (1.33) 

The definition of the n-point spectral correlation implies that they are all symmet- 
ric. Thus, the interchange k' ^ k shows that the above equation is antisymmetric. 
Therefore, we can conclude that equation (1.29) holds. 

As we said, the advective term redistributes energy transferring it from one wavenum- 
ber to another. It has no infiuence over the total energy: 

J 771 poo 

— + / 2iyk'^dkE{k,t) = 0, (1.34) 
dt Jo 

the rate of decay of the total energy per unit mass is the dissipation rate e — dE/dt. 

Henceforth, the advective non-linear term represents the collective action of all the 
modes over a specific one. Its general expression (1.28) can be rewritten in an integral 
form; we use the property (1.33) to put the spectral (2-point correlation) equation as 
follows 

J POO 

—E{k,t)= dk' S{k,k',\k-k'\,t)-2iyeE{k,t), (1.35) 
dt Jq 

where S satisfies the equation: 



pco pk2 

/ dk' dk S{k,k',\k-k'\,t)^0, (1.36) 

Jki Jo 

for arbitrary ki and k2- 

Previously we introduced the Kolmogorov hypothesis for the configuration space 
assuming the process is also time-stationary, but equation (1.35) does not posses that 
property. Let us consider for a moment that the advection is absent, then we have 

E{k, t) = E{k, to) exp [-2uk'^ {t -to)] : 



the greater the wavenumber the faster the energy density will decay. 

So, that is how the cascade happens: the non-linear term takes energy from the 
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low wavenumbers, where there is net energy production, to compensate the net losses 
due to viscosity dissipation at high wavenumbers. We hope that this transfer will lead 
the system at large times to a steady state. Because this situation is found in many 
real flows, the fully developed turbulence model is a representative class of turbulent 
phenomena. 

To consider a time-stationary state within this model we will introduce an artificial 
term. We will restore an external random force- like term f(k.,t) into the spectral 
equation. It should also satisfy, remember equation (1.16), 

k-f(k,i)=0. (1.37) 

It modifies equation (1.27) as follows: 

(J-^ + 2uk'^^ E{k, t) = T{k, t) + 27rk^{i{k, t) ■ u(-k, t)). (1.38) 

Now, in order to find an explicit form for the rightmost term above we must char- 
acterize the random force. Lately, we have argued that the turbulent state should not 
be modified by external forces; then, suppose this force is Gaussian distributed, and 
its autocorrelation is given by 

(f (k, t) ® f (-k, t')) = D(k) W{k)S{t - t'). (1.39) 

While the operator D is introduced to obtain an homogeneous, isotropic, and station- 
ary force, the ^-function makes it highly uncorrelated in time. Finally, W has to be 
described: we will assume that the system response, for small time intervals |t — is 
given by the Green function g{k, t — t'), such that 

u(k,t) = [ dt' g{k,t-t')i{k,t'). 
Jr 

This kernel function also has the property 



g{k,t-t') 



for i < t' 

1 for i = t'; 
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so, it is causal and recovers the acceleration at equal times. We write then 

27rA;2(f (k, t) ■ u(-k, t)) = 2nk'^ I dt' g{k, t - t') tr (f (k, t) ® f (-k, t')) = ink'^W{k). 

Jr 

(1.40) 

Henceforth, the equation (1.38) achieves its final form 

^E{k, t) = T{k, t) + Ank^Wik) - 2ueE{k, t). (1.41) 

Stationary in time is found when the right-hand side of the latter equation is zero, 
under this circumstances it yields: 

poo 

/ dk' S{k,k',\k-k'\,t)+ink'^W{k) -2iyeE{k,t) ^0. (1.42) 
Jo 

If we integrate this equation over the whole k-space we obtain, 

/•oo poo 

/ Ank'^dk W{k) = / 2uk'^dk E{k) = -£, 
Jo Jo 

but a well-posed problem with separated input and dissipation ranges, it is what we 
have in the inertial range, implies the existence of a wave number k* such that the 
former is replaced by 

pk* poo 

/ Ank^dkWik)^ / 2uk^dkE{k)^-e. (1.43) 

Jo Jk* 

This means that the input term is peaked around A; = 0, and that the Reynolds number 
should not be too low. 

Two energy-balance equations can be now drafted from (1.42): for the first we 
integrate from zero to k* 



pk* poo pk* 

/ / dk' dkS{k,k', \k-k'\)+ A7rk^dkW{k)^0, (1.44) 

Jo Jk* Jo 

here we have used property (1.36) to set the integration interval of the advective term; 
the second equation is obtained with the same argument but integrating from k* to 
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infinity: 

poo nk* poo 

/ / dk'dkS{k,k\\k-k'\)- 2iyk'^dk E{k) = 0. (1.45) 

Jk* Jo Jk* 

While the first equation tells us that the energy injected to the system from the 
low modes are transferred by the non-linear term to the higher modes, i.e., the inertial 
forces transfer energy from low to high wavenumbers. The second equation explains 
that the energy transferred is dissipated in the range k* < k' < oo. 

We have finished characterizing the stochastic properties of the turbulence; also, we 
provided a model for the energy cascade. The inertial range is thus defined by those 
wavenumbers lesser than some k*, where most of the advective term is concentrated. 
Therefore from the dimensional arguments we have used in Section 1.1.1, we can take 
kd — k* 27: / Iq We mentioned before that the injected energy W should 

be concentrated around A; = 0; so, the hmit 27r/L is the cut-off that sets the injection 
of energy. This is how the inertial range (see Figure 1.1) is set within the wavenumber 
space: 

271" 

— <:k<^kd. (1.46) 



Now, the spectral density of energy becomes independent of the viscosity in this 
range as long as the Reynolds number is high. There is no other possibility for this 
function than to be 

E{k) = a£2/3fc-5/^ (1.47) 

with a an adimensional constant. But, of course, this is the case of — > which is an 
idealization. If we introduce the dissipation via kd the distribution should be written, 

E{k) = ae^/'k-'/'F (^^^ , (1.48) 

with F(0) = 1. 

The theoretical form for F is still under discussion. Moreover, there are strong 
clues that suggest that this functional relation extends all over the inertial range. It 
will be discussed next. 
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Figure 1.1: Energy balance: three intervals are marked: below 27t/L (in blue) the 
energy is injected by the external forces, above k* (in red) the energy is dissipated 
by viscous forces, and finally the inertial range comes from the balance between the 
energy injection and dissipation — as it is shown in equation (1-42). 



1.1.3 The problem of the intermittency: Kolmogorov refined 
hypotheses 

The phenomenon of intermittent turbulence was first experimentally noticed by Batch- 
elor and Townsend (1949). They found that the energy was nonuniformly distributed 
throughout space in a fully developed turbulence. Some regions showed to be more 
active than others. The energy intermittence also implies that the dissipation behaves 
in the same way: this contradicts the assumptions that lead to equation (1.7). That 
is, the volume VJ occupied by eddies of size / is not constant, as we supposed after 
(1.6). Therefore, the global average should be replaced by local averages of the energy 
dissipation rate, for the former does not represent the behavior of e. 

citetpaper:oboukhov was the first to tackle this problem. His proposal was to divide 
the spatial domain into a collection of ensembles with characteristic dissipation Sr — 
where Er is the locally averaged dissipation over a spherical volume with diameter r 
and center r'. Later Kolmogorov (1962) used this proposal to rebuild his hypotheses. 
He added another hypothesis to shape the randomness of the energy dissipation rate: 
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Sr{r', t) is a Gaussian variable having as variance of log£r the following 

cT^(r',t) = A(r',t) + 9/xlog 

is a universal constant. Eventually he found the following scaling law for the n-order 
structure function: 

(T \ Am(n-3)/3 
-) ' (1-49) 

where Aur = u(r + r') — u(r'), is an adimensional constant, and e is, as before, the 
energy dissipation rate averaged over the entire bath of characteristic size L. Never- 
theless, it is found that the universal constant /i depends on the order n. This is known 
as scale-similar theory of intermittency. It is worth nothing that it is unnecessary the 
log- normality to explain (1.49). 

To get a picture of how intermittence is created from the continous process of 
stretching and twisting of the advective term we will follow Prisch et al. (1978). Let us 
restart from equation (1.5), assume that the average number of offspring of any eddy 
is N — an eddy of scale In is supposed to give rise to N eddies of scale In+i, irrespective 
of the value of n. These N eddies occupy a fraction /3 of their predecessor . This 
fractional reduction in volume from one generation to the next is given by 

P = = Nn^ < 1. (1.50) 

Furthermore, let us suppose that the largest eddies fill all the space available to them 
Vl L^, the n-generation occupies just its active fraction 

Vi„ = P'^Vl. (1.51) 

With the arguments we have given, the accumulated energy from the eddies of size n 
is now 

Ei^ - Vi^vl - VlP^'vI (1.52) 
while the globally averaged energy flux gives — remember that £ ~ because there is 
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no energy losses within the inertial range — the following scaling law for the velocity 

£ 77 77 P Vi^. (1.53) 

However, we want to express the above formulas in terms of the scale length. Let 
us assume that the splitting number N is inversely proportional to a power of the 
contraction ratio, i.e. N — n"^. The volume reduction for the n-generation of eddies 
results ^ 

(3- = (tVtt^)" = ttM" = in"f-^ = (^j-J , (1.54) 

here we have used definitions (1.50) and (1.5). Now, we can estimate the intermittent 
scaling law for the energy 

3-D 

EJV, ^ e-lHT (I) ' , (1.55) 

D is identified as the Hausdorff fractal dimension by Mandelbrot (1974), and represents 
how much space is filled by the eddies so 3 — D > 0. The functional expression for 
2-order structure function, or just structure function, should be 

(||Aur|n = X„(sr)^/^(9^, (1.56) 

and comparing against equation (1.49) it is 2// = 3 — 

Intermittency means that the probability for having large velocity fluctuations in- 
creases at small scale, vi^ ~ s^l'^tn^iLjl^^^'^^l'^ ^ but at the same time the amount 
of these eddies decreases with the scale. Furthermore, this imphes the existence of 
a singularity in the Navier-Stokes equation when i/ — > 0. It can be stated formally 
as follows: given h — 1/3, these singularities are contained in a set Sh with fractal 
dimension Dh — D. In the sense, r' e Sh 

lim^^O. (1.57) 



Therefore, we have proven that the Frisch's /3-model accounts for all these prop- 
erties of the intermittence, and produces a 2-point structure function coherent with 



1.1 Turbulence 



19 



the experimental findings. Besides, we can estimate higher structure functions for the 
velocity field, 

dlAurll'^) ocr^" (1.58) 

with (n = hn + 3 — D, and h = {D — 2)/3. The structure function exponent has a linear 
growing with n; actually, experimental tests (Anselmet et al., 1984) shows a non-linear 
grow. The statement affirming that the energy transfer per unit volume is constant 
through the different eddies scales is not valid (Frisch and Parisi, 1983). 

To override this problem let us suppose that instead a single scale h there is a 
range of them hj^in < h < hmax where (1.56) happens. Thus, we have a wider set 
of singularities, ^hG[h^in,hmax]^h- Each subset has its own fractal dimension Dh. To 
calculate the moments of the velocity variation Aur at a point r' we have to look for 
the probability of finding such variation in one of the subsets Sh- This probability is 
proportional to the volume of Sh thickened along the normal direction by a length r. 

This enlarged set is defined for any given set F C M" as the A-parallel body: 

Fa = {r e M" : ||r' - r|| < A, r G F}. 

We want to estimate its volume, let us pick up some examples first. For a single point 
set, F = {r} it is obviously voI(Fa) = 47rA^/3. If we take F as a segment of length / 
it is vol (Fa) ~ tt/A, and for an extended flat surface of area A is vol (Fa) ~ AX^. In 
each case, voI(Fa) ~ cX^'^ with s the dimension of F. This idea can be extended to 
fractal dimensions (Falconer, 1990). It follows that the n-order structure function is 
the average 

(IIAurD oc y"/i(rf/i)r"'^+3-^'' (1.59) 

where ^{dh)r^~^^ is the probability over the spectrum [/imin, ^max]- In the short dis- 
tance limit this integral can be estimated using a saddle point approximation, so 

(II Aurir) with Cn = min(n/i + 3 - F>,,). (1.60) 

h 

Since intermittency is related to singular velocity variations one should expect h < 
hraax = 1/3. Of course < Dh < 3 and then ^„ < n/3. Also, we can see that each 
depends on a specific value of h; therefore, the velocity moments pick out a particular 



1.2 Passive Scaleir Fields' Chciracterization 



20 



subset Sh- 

Besides the latter model, other alternatives to the /3-model using fractal geometry 
has been proposed (Benzi et al., 1984) to explain the energy transfer: the random 
/3-model. It is assumed that the eddy splitting is not constant; from one generation 
to the other N changes. That is, the contraction factors /3 are independent random 
variables. 

These fractal models gives a better understanding of the phenomenon of advection 
and intermittence. But neither of them predict the values the anomalous dimension, 
and it must be found by experimental means. 

1.2 Passive Scalar Fields' Characterization 
1.2.1 Scalar turbulence 

The turbulent flow transports and disperses any scalar by making parcels of fluid 
follow chaotic trajectories: the non-uniformity of the turbulence causes lines of constant 
scalar stretch and fold. This process drives the scalar concentration through smaller 
scales; eventually, the diffusivity k associated to the scalar (e.g. thermal diffusivity, 
molecular diffusivity, etc.) prevails over the advective mixing. There are so many 
parallels between this behavior and the one of velocity field that we address it as scalar 
turbulence (Shraiman and Siggia, 2000). 

In a turbulent flow the scalar is controlled by two processes: transport, the physical 
translocation of the scalar via the combined action of fluid advection and diffusion; 
and mixing, the irreversible decay of fluctuations because of the scalar diffusion, k, 
that tends to reduce the scalar field to uniformity. In a turbulent flow, both processes 
become independent of k as it goes to zero. This limit is the so called fully developed 
scalar turbulence. 

Any given scalar (concentration of the) quantity © put into a static fluid is subject 
to a diffusion equation, that is, 

^e(r,i) = /.Ae(r,i). 
The extension to flowing fluids is accomplished replacing the partial time derivative by 
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the total derivative, so to speak: 

e(r, t)^^ e(r, t) + {v V)e(r, t)^K Ae(r, t), (1.61) 

with V the random velocity field which is naturally solenoidal, i. e.,V • v = 0. 

Now, the above equation is our starting point. One should consider introduce the 
Navier-Stokes equation to describe the velocity field, but we will not do so. In fact, most 
of the actual developments in scalar turbulence does not need a deep understanding 
of the velocity field: for all the purposes here will be enough to describe it as a given 
isotropic and homogeneous stochastic field. 

Because of the incompressible nature of the fiuids equation (1.61) can be rewritten: 

4-© + V-(ue-«;Ve) = 0, (1.62) 
at 

here u is the stochastic velocity field for the zero mean turbulent velocity field, and 
we also used V • (u ©) = u ■ V0 + (V ■ u)0 — with the last term vanished. Wc write, 
again, the concentration of the scalar as the sum of a mean and a random fiuctuation, 

e(r,t) = (e)(r,t)+i?(r,t), (1.63) 

with = 0. In this case the mean of the scalar (©) can not be neglected, because 
scalar turbulence is usually generated by maintaining a mean scalar gradient. Note that 
the homogeneous and isotropic velocity field does not guarantee the same properties 
on the scalar turbulence. 

Henceforth, using the latter definition, the averaged equation (1.62) gives us, 

4- (0) + V-((u^?)-KV(e)) = O. (1.64) 
ot 

The first term within the divergence term represents the effect of velocity field in the 
transport of the concentration. We may assume that the bulk effect of this transport 
is proportional to the action of the mean concentration's gradient. That is. 



(u^?) = -«TV(e), 



(1.65) 
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kt is defined as the eddy diffusivity. This prescription for (wd) is widely used in 
engineering apphcations, for it 'solves' the closure problem. It also expresses that 
anisotropics introduced by the large scales are maintained overall the scales range. 
Henceforth, the universality of the Kolmogorov treatment could be lost, because the 
turbulence can not be unbounded from the large scale forces. We will discuss the full 
effect of this behavior in forthcoming sections. 

The integral 1/2 d^r ■j?^ is a good measure of the concentration of the scalar field. 
It is zero only when i? is zero, and therefore no turbulence is present. An energy balance 
equation can be obtained for it: we just take the difference between eqs. (1.62) and 
(1.64), and afterwards the average to have. 



integrating it over the whole bath, where the divergence terms are suppose to vanish. 



We have used definition (1.65) to arrive to this last equation. The scalar will be 
stationary when d{d'^)/dt = 0. Then, there is balance between the injected energy 
from the scalar gradient and diffusive term: 



while the left-hand term represents the scalar concentration created per unit time, the 
right-hand is the concentration destroyed by molecular diffusion. 

But if we look for isotropic and homogeneous scalar turbulence the gradient of the 
mean field must be zero, and also V ■ (wd) =0; therefore, the former model for energy 
injection should be left behind. Time-stationary turbulence needs an energy source to 
counterbalance the scalar dissipation, so an external 'force' / should be supplied at 
the right-side of (1.62). Doing so, in equation (1.67) we just make the replacement 



yields 




(1.67) 
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Let us inspect briefly the decaying of the scalar concentration. It is appropriate to 
use the Fourier formalism we applied with the velocity fleld; from equation (1.62) with 
the above prescriptions 

^-ka) ^{r,t) = -V • (ui?)(r,i), 



9t 

we get the following representation: 

(4- + ^^^] ^(k, t)^-i f d^k' k • u(k - k', t) i?(k', t). (1.69) 

Obtaining an evolution equation for the non-stationary spectrum is straightforward: 
multiply both sides of the above equation by i?(— k, t), and then average. The result is 
the following balance equation for the spectral distribution of the concentration: 

{j^ + Kk^^F{k,t)^T^{k,t), (1.70) 

where it is defined F{k, t) = ^nk'^ t) T?(-k, t)) with 

/■oo 

/ dkF{k,t)^{'&^), (1.71) 
Jo 

and the scalar transfer spectrum 

T^{k,t) ^ -87rk^ I d^k' ik- {^{\i-\s!,t)n{\s!,t)^{-\i,t)). 

Now, as in Section 1.1.2 it can be interpreted in the same way as the energy-balance 
equation. The scalar transfer also possesses a conservative property: 

/•oo 

/ dkT^{k,t)^0, 
Jo 

this is because of the isotropy we have imposed on the transport term (wd). Taking 
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integrals on both sides of the spectral balance equation we obtain again, 

^ Jt^^'^ = -X = I k'dk F{k,t). (1.72) 

The idea that all the properties encountered studying the turbulent velocity can be 
mapped to the scalar turbulence was also employed to build an inertial convective (time- 
stationary) spectrum for the scalar variance. Oboukhov (1949) and, independently, 
Corrsin (1951) settled down the first steps towards a scaling law for scalar fields. Using 
the analogy between the scalar diffusion coefficient k and the viscosity coefficient they 
both used dimensional arguments to find an extra cut-off for the energy spectrum: 

A;, = £V4^-3/4. (1.73) 

Because we now have two ultraviolet cut-offs, and k^, the inertial convective 
range should be that where both the viscosity and the diffusion coefficients go to zero, 
i.e., 27r/L <^ /c -C mm {kc, kd}- Afterwards, within this range the power spectrum for 
the scalar is, 

F{k) = Axe-^^^k-^/^, (1.74) 

where A is known as the Oboukhov-Corrsin constant. 

Unlike the velocity field, there is an innate term (u •&) injecting energy down to the 
small scales, and it can not be discarded at will. It causes one of the characteristic 
features of the scalar turbulence: the observation of ramp-and-cliff structures regardless 
the model introduced for the velocity field. This is the source of scalar intermittence 
and anisotropy seen in experiments and numerical simulations. On the contrary, the 
Kolmogorov-Oboukhov-Corrsin (KOC) theory assumes that the advective term restores 
universality, i.e., independence of large-scale injection mechanisms, and thus isotropy 
at the inertial-convective range. This is far from our final goal, which describes the 
full behavior of the scalars. But before inspect the anomalous scaling inherent to 
intermittence, we will introduce an isotropic model binding the scalar fields to the 
velocity scaling behavior. 
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1.2.2 Kraichnan's model 

The Kraichnan's model (1968) for passive advection assumes that both velocity and 
scalar fields are homogeneous and isotropic. So, we must introduce an external force 
/ — proportional to the diffusion k, — in order to compensate the energy dissipation. 
Equation (1.62) is modified to yield, 

^ e + V- (uG-Av Ve) =/. (1.75) 

Let the velocity field be a Gaussian process independent of the random force. This 
synthetic velocity has the following 2-point correlation 

(u(r, t) ® u(r', t')) = D(r - r') S{t - t') (1.76) 

with D(r) = D(0)I - d(r) such that. 



d(r) = D 



IrIP 



r||<, (1.77) 



where < C < 2 is a fixed parameter, and d{— 3) is the dimension of the configuration 
space. More rigorous approaches require a regular covariance function, so it is often 
introduced an infrared cut-off m to rewrite it as follows 

where is related to the constant in (1.77) by 

D = r[(2-c)/2] 

22+C7r3/2^(3 + C)r[(C + 3)/2]''°' ^'"'^^ 

as a direct calculation shows. 

This model for the stochastic velocity field is far from realistic. Besides the fact 
that it mimics the spatial dependence of the correlation function, the opposite happens 
with the time. The velocities are white noise in time: they are uncorrelated everywhere 
but t — t' . Hence, at any instant all the moments of the velocities are infinite! There is 
independence from the past, which contradicts the KOC model where the time-to-live 
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of the scalar inhomogeneity is T(r) ~ r^/^. 

Moreover, the KOC theory predicts the (||Au|p) ~ r^/^ ~ r^/^/T(r) law. When 
compared against the Kraichnan's model (||Au|p) ~ r'^S{t): we have C = 4/3. It is 
said then that the velocity field changes very rapidly in time. Additionally, we observe 

that e»(||u||) ~ VF^. 

Also, we will set the force to be Gaussian with zero mean, and having the following 
2-point correlation function: 

{f(v,t)f(r',t'))^C(r-r')5(t-t'), (1.80) 

where C is invariant under rotations and its compact support set has the extension 
of the bath L. We are interested in the stationary regime which is found when (1.68) 
holds: it requires 

X^K{{V^r)^{f^)^lc{0). (1.81) 

Now, with these prescriptions we will solve equation (1.75). Prom the homogeneous 
equation, which is a Fokker-Planck equation, we obtain the Green's function: 

G,{v,t;v',t')=0 (1.82) 

with initial condition Gk(i", ^o; i"', ^o) = ^^(r — i"'). Hence, the Kraichnan's equation has 
the solution, 

^{r,t)= f dt' [ d\'G,{r,t-r',t')f{r',t')+ [ d\' G,{r,t;r' ,to)M^'), (1.83) 

Jto JR^ JR^ 

for^(r,to) =^9o(r). 

The resolvent is directly found from the Lagrangian trajectories of the fluid parti- 
cles, that is, 

dp{t) =u{p{t),t)dt + \/2KdB{t), with x (to) = xq 

where B the isotropic Brownian motion — note that both terms are of order \/t. This is 
the symbolic representation for a Stochastic Integral, and we will delay its interpretation 
until Chapter 3. Instead, we are going to follow the original Kraichnan's paper. There, 



d_ 
di 



+ u • V - K A 
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the resolvent is found to have the formal expression, 

^ ' . (1.84) 



^^(r, t; r', to) = exp | - ^ dt' [kA-u- V] {t') | 



(r,r') 



If we want to estimate the solution at time tQ + St, we just need to approximate up to 
the first order in St of the Green's function, that is. 



G^{r,to + St-r',to) = 1+ / dt' {u{t')-V\^y^ - KASt\^^y^ + 

J to 

Jrto+St ft' 
' dt' dt"(u(O-V)(u(O-V)|(r,r0 + 

to to 



This procedure allows us to build a differential equation for the moments of the scalar. 
But here, the first order approximation let us only estimate the time evolution of the 
2-point correlation: 

(i?(r, to + St)i»(r', to + St')) - (^?(r, to)i»{r', to)) = « ( Ar + A^O (^?(r, io)^(r', ^o)) St+ 
rto+St ft' 



+ 



rto+ot pt' 

/ dt' / dt"{{n{v, t') ■ V)(u(r, t") ■ V) ^(r, to)^{v' , to)) + {r ^ r'} + 

Jto Jto 

pto+5t i-to+St 

+ dt' cii"((u(r,0 • V)(u(r',0 • V)^(r,io)^(r',io))+ 

Jto Jto 

rto+St pto+St 

+ dt' dt"{f{r,t')f{r',t")). 

Jto Jto 



The process t?(r,to) is markovian^, for it is statistically independent of u{t) and f{t) 
whenever t > to- Therefore, assuming the scalar is time-stationary the left-hand side 
is zero. Under the homogeneous hypothesis it is F2{r) = {'d{r)'d{0)) and so we have: 

(-K A - V • d(r) • V) F2(r) = C(r), (1.85) 

adding the isotropy condition implies the operator — ^ " d(r) • V can be rewritten 

^It is said that a process Xt is markovian or possess the Markov property if the future behavior 
of it given what has happened up to time t is the same as the behavior obtained when starting the 
process at Xt (a detailed description can be found at Shiryayev (1984, Chapter 7)). 
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as 




Once we set the boundary conditions, the solution to (1.85) is: 




.d-l' 

1 



With a bath having finite extension L, the inviscid hmit k, ^ yields 



+ ..., 



(1.86) 



Dd{d -1){2-C) 



the constant C comes from the selected force correlation C(r). Besides, the 2-point 
structure functions is independent of it, and is 



Therefore, homogeneous and isotropic scalar turbulence imposes a universal law for 
the correlation function, since it only depends on the mean dissipation rate and the 
distance. 

1.2.3 Anomalous scaling, anisotropy and diffusion 

We surveyed with the Kraichnan's model the behavior of the isotropic and homogeneous 
scalar turbulence. We note that intermittence in the velocity field is directly translated 
to the scalar field, equation (1.87). The anomalous scahng in scalar turbulence is just 
the separation from this inherited power spectra. 

In the present context universality of scalar fields should be understood as the 
existence of limiting correlation functions (n„i?(r„, i„)) in a stationary regime when 
K,m,l/L — > — as defined in the former section — are independent of the external 
force, i.e., the shape of the correlation function C. This condition let us build the 
following inertial range: 



4x 



(1.87) 



Dd{d -1){2-C) 
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with D as in (1.77). 

Gaw§dzki and Kupiainen (1995) determined that the correlators become indepen- 
dent of the diffusion and infrared cut-off, but there exits a dependence upon the external 
force. The integral scale L contributes to the n-point structure functions 

S2n{r) = ([^(r) - ^(0)]^") ~ A||r||^" = (L/||r||)''- ||r||(^-^)^ (1.88) 

for ||r|| <^ L as K,m. go to zero. The amplitudes 72n depend on ( and the function 
C, while the exponents p2n are just functions of C but not C. While for n = 1 the 
intermittence is absent, p2 = 0, for n > 1 the p's are positive, increasing, and convex 
function of n. 

They followed the arguments from Kraichnan (1994) who noticed that an exact 
expression can be written for advection effects on scalar structure functions when ve- 
locity fields change rapidly in time, they are delta correlated. The second assumption 
is the need for the scalar field to be Gaussian, because it makes possible to break 
down the closure problem, that is, the dependence on higher moments of the scalar. 

Next, the velocity statistics is introduced through the two-particle eddy diffusivity 

V{r)^l: I rft'((5p(r,i)(5p(r,0) 



2 „ 

where S\\u{r, t) — [u(r', t) — u(r' -|- r, t)] ■ r/r. Inside the inertial range it has the form 

v{r)^Vo{^) , (1-89) 

here we have changed the former exponent ( into a functional of the diffusivity C(^)- 
The final differential equation for the n-point structure functions is 

dS2n{r) 1 d f , , dS2n{r) \ _ Ac^ m^'^^-^^) 



dt r'^ dr \ dr J S2{r) 

The steady state makes the temporal derivative vanish, and the 2-point structure func- 
tion is easily obtained: 

^ ^^'^2(0) / /^^C2W n 90) 
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where kAS2{0) is the rate of dissipation of the scalar x, and the scalar exponent is 
related to the eddy diffusivity exponent by 

C2(^) = 2 - Civ), with < Civ) < 2. (1.91) 

These equations guarantee a precise functional relation for the exponents within the 
inertial-range, and again this exponents relation is found independent of the external 
source: 

C2nW = ^Vl2<2W + [3-C2W]^ - ^[3 - C2 W] (1.92) 

providing an asymptotic behavior as n increases C2n ~ \/3n^. If we suppose that the 
above equation can continuously extended to n = 1/2 we found then the exponent of 

{\A^\) to be 

Ci(^) = IVocMW^U&W- ^[3 - C2(^)]. 

Constantin et al. (1991) showed this exponent is also related to the fractal dimension^ 
of the set of isoscalar surfaces contained in a sphere of radius the order an inertial-range 
scale, that is, 

dimH^'\c) = 3 - Ci(i?). (1.93) 

It resembles our discussion of the volume fraction P filled by eddies with characteris- 
tic length In in Section 1.1.3, but in this case the scalar exponent is exactly the co- 
dimension. Meaning the scalar field exponents are strongly determined by the scalar 
spatial distribution. 

Furthermore, we can examine the significance of the exponents in the limiting cases. 
As Civ) ~^ 2, all the structure function's powers go to zero as (1.91) and (1.92) shows. 
This is reahzed when the scales are near the inner length defined by the velocity field, 
where the viscous term in the Navier-Stokes equation becomes relevant. We note that 
dimj^ ^9~^(c) — >■ 3 and the scalar fills all the space. 

The opposite limit, Civ) ~^ 0) makes the power spectra approach to k~^. But at 
the same time the effective eddie diffusivity rio grows. Thus, we observe from both 
equations (1.90) and (1.87) — with D from (1.79) — their coefficient going to zero. 

Because at small r the difference | Aru| is at most of order 1, it is < C2n < 2n. Also, 

^The fractal dimension dimn is formally known as Hausdorff dimension. See the Appendix A for 
a definition and some relevant properties. 
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applying the Holder inequality we can prove Q/q < Cp/p, for any two p,q > 0. So we 
found an upper bound to the scalar dimension, as defined in (1.93), 2 < dimff '&~^{c) < 
2 + ({r])/2. Henceforth, in the limit ({r]) — the scalar turbulence is contained in 
2-dimensional sheets. 

Additionally, if the velocity field u changes slowly in time, but remains Gaussian, 
and has a long scahng range < C2(u) < 2 — in the sense of (1.58) — it acts hke a rapidly 
changing field because the large scales sweep fiuid elements rapidly through the small 
scales. So we should have C(?7) = C2(u) + 1, but it does not hold for any pair of values 
(■(77), C2(u). As can be seen from the ranges covered by each power; moreover, it seems 
valid near the classical values: C,{rj) — 4/3 and C2(u) = 1/3. In any case we can affirm 
that there is not exist an invective relation between both variables. 

Up to now, we have seen that the effect of the external source can not be separated 
from the phenomenon of intermittence. Besides, scalar turbulence is generated from 
the existence of a mean field of the scalar. So the prescription of an isotropic external 
source must be abandoned: the scalar fiuctuations are excited by entwining of the mean 
external gradient of the passive scalar by turbulent fiow. 

Therefore, given the external gradient V(©)o 7^ equation (1.62) yields. 



The transference of the scalar concentration is done through the wavenumber space. 
So, let us consider the spectrum of the anisotropic scalar in the wavenumber interval 
27r/L < k < mm{kc,kd}- Because the molecular diffusivity is supposed to go to 
infinity, the two relevant terms are the advection and mean scalar gradient in the 
above equation. Both must be of the same order. 

Besides, the eddy diffusivity, as we briefiy introduce it before, measures the rate 
of transport of scalar concentration from one portion of fiuid to another. It was first 
introduced by Heisenberg (1948) (see McComb, 1991, p. 75) and as a function of the 
wavenumber is written. 



If we suppose the energy is peaked around A; = 0, the most important contribution to 



d_ 
di 



v.(u(e)o). 



(1.94) 




the diffusivity is around k: thus, r]{k) ~ ■\jE{k)/k 




1/2 
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is the magnitude of the stochastic velocity field at a given wavenumber for 

/•oo 

(||u|p),= / E{k')dk'. 

J h 

Henceforth, given the eddy diffusivity we can replace the advective term by r]{k)A'& 
and obtain 

u- V(e)o ~ r](k)A'&. 
So we estimate the scalar spectrum with this equation and using definition (1.71), 

F{k) ~ fc-^||V(e)of . (1.95) 

It was encountered that a renormalization procedure (Elperin et al., 1996) provides 
this spectrum for the anisotropic source. Also, the isotropic case is undertaken with 
this formalism. Given the power law ^j^g velocity spectra it is found 

F{k) oc A;-2+^2(")/2 (i gg) 

whenever the quotient between the effective viscosity u{k) is greater than the kinematic 
uq. In particular this is condition is broken when C2{u) 0: it is i/(A;) \ uq. 

We are seeing two processes in this discussion: the largest eddies grabbing energy 
from the anisotropic mean field, they significantly contribute to the structure function 
below certain scale l/A;*; above this characteristic length the external forces turns 
isotropic introducing the suitable power spectra. 

These anomalous exponents sensibly modify the probability distribution of the 
scalar fields. That is. they deviate from the Gaussian shape adding (exponential) tails. 
But because in this work we are just interested in the second moments the Gaussian 
approximation is enough. 

Among all the passive scalars fields, the kind represented by the temperature re- 
quires more attention. Usually its turbulent state is reached through convection, but 
because temperature differences are the trigger for turbulent mixing the temperature 
field should not, at first, be passive. Whether the temperature is active or passive de- 
pends on the magnitude of the contribution from buoyancy forces to the total energy. 
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It is despicable when the Rayleigh number is small: 



7^a = 



agATL^ 



«1, 



(1.97) 



where a is the volume expansion coefficient, g is the acceleration due to gravity, and AT 
is the temperature difference between bottom and top of the bath of size L. Meaning it 
does not affect the state of turbulence. It happens when this number becomes relevant 
that scaling and gaussianity of the temperature field are observed (Sano et al., 1989; 
GoUub et al, 1991; Ching, 2000) (for the soft convective turbulence TZa < 10^). But 
the scaling exponent behavior is more complex than in the scalar turbulence problem, 
and most of the arguments presented here are not applicable. 

1.3 The Turbulent Index of Refraction 

We are going to finish this chapter describing the behavior of the refractive index inside 
a turbulent fiow. It has been proven long ago (Bean and Dutton, 1968) that 



where it is: n the index of refraction, T the absolute temperature, P the air pressure, 
and e the water vapor partial pressure (both in millibars). This equation is consid- 
ered valid for frequencies ranging from IMHz to around 30GHz or more. For light 
propagation it is assumed that the humidity term is negligible, hence 



Without buoyancy effects present we can assume the local pressure to be constant; 
thus. 



That is, deviations from the mean temperature field are proportional to those of the 
refractive index. Adiabatic corrections must be introduced to this formula when the 
system is the whole atmosphere (see Ishimaru, 1997, p. 523), but they do not break 





(1.98) 
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the hnear relation. Therefore, the deviations from the mean fields are proportional, 
and the refractive index inherits the passive scalar properties from the temperature 
field proved it is a scalar field. 

As we introduced a synthetic turbulent velocity field to make insight into the be- 
havior of the scalar turbulence. We will do the same for the lightwave propagation in 
turbulent media. That is, we will introduce a model for the refractive index such that 
most of the properties described in the former sections are present in it. 

Our assignment is to associate to the turbulent index of refraction a suited stochastic 
process. We have seen that the anomalies in the exponents of the n-point correlation 
functions drives the scalar fields apart from the Gaussian statistics at first. Never- 
theless, all the theoretical models presented here induce a Gaussian behavior for the 
scalars because the stochastic velocity field is Gaussian. Moreover, because in this work 
we will only be interested in the first moments of the propagated light, the Gaussian 
distribution is enough. 

The family of Gaussian processes is wide, and each member of it is defined, as 
we shall see in the next section, by its 2-point correlation function or covariance. 
Because we are specially interested in the properties of the atmospheric propagation, a 
fully developed turbulence, the covariance defined within the inertial range characterize 
almost all its properties. Although until now we have not explicitly given an expression 
for this covariance, we will show that the structure function is sufficient. 

When the velocity field is homogeneous and isotropic, and the external source is 
isotropic the structure function for the stochastic refractive index, according to (1.90), 
should have the form: 

Mr + r') - //(r')]') = A^M^, k < ||r|| < L, (1.99) 

/J, is the stochastic component of the turbulent refractive index n, the constant A2 is 

known as structure constant^, and < C < 2 . Even the limit C ^ 2 must be taken in 

consideration because of the anisotropic behavior of the passive scalars 

The structure function is in fact the covariance of the increments. The translation 

invariance of equation (1.99) implies that these increments are stationary in the statis- 

^Usually is noted as when ( = 2/3. Its range is around 10~'^'*-10~'^"'^m~^/'^ for low altitude 
measures and 10-i*-10-^'^m-2/3 for the high altitude (Tatarski, 1961; Ishimaru, 1997). 
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tical sense. That is, the probabihty distribution remains the same under translations. 

On the other hand, previous works make different prescriptions for the turbulent 
index from the ones we have given. They assume the process should be stationary and 
give other covariance functions in consequence, for example: the exponential correlation 
function, 

{/.{v + v'Uv')) = (/.2)exp(-||r||V/^) , (1.100) 

used by Bcckman (1965) and Ishimaru (see 1997, p. 358), to solve some propagation 
problems; or the uncoupling relation 

{l^{z + z',p + p')i^{z', p')) = 5{z)A{p), (1.101) 

that makes the process markovian in the propagation direction z, first suggested^ by 
Klyatskin and Tatarski (1970) and later mathematically formahzed by Leland (1989). 

But we have seen the stationary property for the stochastic refractive index is 
unnecessary; moreover, the 2-point correlation function (1.99) only demands stationary 
increments. Our proposal will be to use the fractional Brownian motion as a model 
for passive scalar fields (Perez and Garavagha, 2001). 

1.3.1 The synthetic refractive index: the fractional Brownian 
motion 

Finally, in this section we will construct the synthetic index we are going to employ 
studying lightwave propagation. For such a task the refractive index must follow most of 
the properties described above. From all the possible Gaussian processes, the fractional 
Brownian motions seem the best suited to accomplish this as they present the following 
second moment of the increments 

E[(S^(i) - B"{s)f] ^\t- sf" , (1.102) 

in the 1-dimensional case. On the other hand, the equation (A. 14) provides a simi- 
lar expression for the tridimensional space, E[(i?^(r) — i?^(r'))^] = n^=i ~ ^if^- 
This is anisotropic or nonstationary in statistical terms; therefore, it can not fulfill our 
*see the Appendix B for a short description of this model. 
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prescriptions for the Structure Function. Instead of this n-dimensional version we will 
use the change of variable property (A. 12) to introduce the following Gaussian process 

S»:=S^(||r||). (1.103) 

We will call it isotropic fractional Brownian motion. The variance of its increments is 
given by 



E 



S^(r + r') -^^(O 



|r + r'||-||r'||)'^, 



and when ||r|| ^ \\r'\\ or ||r|| <^ ||r'|| we can approximate this equation by 



E 



2 



S^(r + r') -S^(r') ^ ||rf ^, (1.104) 



and so the increments are locally stationary. 

Now, the departure /i from the mean index of refraction tiq is a small quantity. 
That is why the stochastic index also measures the behavior of the permitivity, i.e., 
£(r) ~ + 2//(r). It is a convention among the hterature to substitute the turbulent 
index by the stochastic permitivity e(r) = 2//(r), and so we will do here. If we write 
the exponent in equation (1.99) as C = 2if, from the former property (1.104) under 
the conditions given, thus we define 

e(r) := aB^{r/l) (1.105) 

with a an adimensional constant and I some characteristic scale. When we consider 
a fully developed turbulence, isotropic and homogeneous, the Kolmogorov hypotheses 
sets H — 1/3. We are considering departures from this ideal situation sol/3<i?<l 
will be our working range. 

Let us compare the structure function of the permitivity against the structure 
function generated by this synthetic permitivity. Thus using equations (1.99), (1.104) 
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and (1.105) we have, 

2 

(1.106) 

we used the seh-similarity property in the third hne. The comparation allows determine 
the coefficient a = 2l^\/A^. Therefore, we must determine the characteristic length 
and the constant A2 of the structure function. There are two physically distinguishable 
scenes for setting the scale /, whether we are in the persistent, 1/3 < H < 1/2, or anti- 
persistent, 1/2 < H < 1^ case. For the latter continuity conditions for the limit H ^ 1 
(Sirovich et al., 1994) set / = /q and A2 = Al^^'"^^ ^ the former just keep / = L and 
A2 = A'L2/3-2^^. 

Furthermore, using the probability density (A. 2) with the conditions given above 
we also note that 

E|e(r) -e(r')| - ||r-r'||^; 

thus, in our model we have C2 = 2i? and Ci = H. 

Although, the relation between these two is linear, and thus does not coincide with 
the Kraichnan's model. We will proof next it is well defined; that is, it is consistent 
with other well known quantities. In particular, the fractal dimension associated to the 
isoscalar surfaces of the isotropic refractive index — contained within a sphere of radius 
Iq: that is, dimife~^(c). 

This isoscalar surface can be expressed as the set 

e-i(c) = |r e : S^(r/io) = a-^c} 
from definition (1.105). If we take the plane set P^-ic = {{^, a~^c) : r e M^} C x M: 



>Sf(r) =4A2||r| 



\2H 



~a2;-2^||r||2^, 



e-'{c) = graph nPa-ic- 
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Moreover, using the properties (A. 17), (A. 18), (A. 19) and considering that 



-3 



dimfiPa-ic = hm = Hm — = 3; 

5^0 -log (5 5^0 - logS 



it is 



dimjj (^graphS^ n Pa-ic) = dimjj graph - 1. (1.107) 



It is necessary, thus, to calculate the Hausdorff dimension associated to the set graph . 
We will accomplish this in the following paragraphs. 

First note that the 1-dimensional fractional Brownian motion is X-Holder continous 
(Falconer, 1990, p. 246), that is, for all < A < ii", 

\B" (t) - {s)\ < M \t - , for some M. 

We observe that our isotropic fractional Brownian motion is also A-Holder continous. 
From the triangle inequality, | — | < Ik ~ I'll) latter equation we thus 

have for any < X < H: 

B"{r)-B"{r') = 5^(||r||) - S^(||r'||) < M| ||r|| - ||r'|| |^ < Mjjr - r'||^ (1.108) 

for some M. 

Let be B^ : [0, 2lo\^ — > R, that is, it will be restricted to a box-set containing the 
sphere Bi^ of radius Iq. From all the coverings to the graph let us choose a box-covering 
{Bg} such that dia,m{Bs) < S — with side of length less than 5/2. It is, 

J^^(graph5^) < J2 (diam(S5))" < NgS', 

where Ng is the number of boxes touched by the graph of B^ . Over the domain we 
have at most loS~^ + 1 boxes. On the other hand, it is 

\B^{S/2, S/2, 5/2) - B"{0) \ < MS^ 

from equation (1.108). Therefore, we have at most MS'^/S + 2 boxes piled at a given 
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box on the domain; so, 



Ns5' < M'5^-^+' + 25'. 



Now, as 5 — * the Hausdorff measure remains bounded if and only if s > 4 — A. So 
from property (A. 16) we have 



dimj^ graph < A — H. 



(1.109) 



Therefore, the fractal dimension can not exceed A — H. Next we will find a lower 
bound to the dimension. We will apply a flavor of the potential theory commented at 
the Appendix A. That is, we will look at the integral 



I E 



\\x — y\\^ 



for the stochastic mass measure ii^j- If we flnd this integral flnite then the fractal 
dimension of the set we are studying will be greater than s. We choose the occupation 
measure: 



it counts how many points of the set B are in the graph of B^. For simplicity, let us 
take the sphere Bi^ centered at the origin. We have 



^^'^{dx)nZ{dy) 



|r-r'||2 + 



B^{v/h)-B^{v'/h) 



-s/2 



d^r d^r' 



Because B^{v) — B^{y') ~ A/'(0, |||r|| — ||r'|||) and the triangle inequality, we have the 
following inequality. 



P|||i 



W < III N _LH|| P||r-r'||(^), 
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for the probability densities 



PlIMI-llr'lllW 



Therefore, 



2 / I II II II /II |2-ff 

exp —z I |||r|| — ||r III 



27r|||r| 



/II |2ii" 



and p||r-r'||(2;) 



exp(- 



-z^ Wx - r 



/||2i?> 



V27r||r-r'||2^^ 



/||2 



r — r 



< (1/2) 



i?^(r//o)-5^(r7/o) 



Hi I 



2 X -s/2 




r — r 



d^r (fr'dz 



r|| _ ||r'|||^(||r-r'||2 + z2)V2 
< al/2) f d\ f d^u f dz puiz) ^7 

J 7m ^ ^ |||u + v||-^;|''(i2K2 + ^2)./2 

Puiz) u^v?du 



0<|M|<1 ||u+v||<l 

1 r-l 



< Z^27r^ / ciz / v^civ / dx 







1 ^0 



m — — t; 



'^(/2m2 + ^2)./2' 



in the last inequality we have used \u — v I < ||v + u||. Since < sJivxY + (1 - v'^) 
vx <2 we change the former inequality to 



r — r 



/||2 



< I^At:^ I dz I v^dv 



2 , A-*/2 



\\u — v\ — {IqU^ + 2;2)V2 
p„(2;) dz \ { H v'^dv 



< 



(1-H) 



(/2^i2 + ^2)V2 

/ u^du [ pu{z){iy + z^'/^dz. 
Jo Jm. 



\\u — v\ — v\ 



2 /•2 



Now, the probability density Pu{z) is bounded in any closed interval of u. Let us 
subdivide [0, 2] into intervals of the form [c"^^5, c"5) with c < 1 and c^^5 > 2 for some 
k, S such that U'^f,[c"'~^^d, c^S) D [0, 2]. We can rewrite the right-hand side of the above 
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inequality as, 



\r-r'f + 



5^(r//o)-5^(r7/o) 



-s/2 



12tt 



2 °o rc"-S 



(1 ~ ~ ^c"+i5 

Making the change of variables u — c^t we have: 



|r-r'||2 + 



B"{r/k)-B"{r'/lo) 
I2n 



2 \ -s/2 



< 



JcS Jm. 



Since the fractional Brownian motion is self-similar we have: 

/ pMrnic'^-t' + z'r^'dz = f Mz){iy-t'+c'-''z'r/'dz 

Jr Jr 

<Mo I {ll^^e + ^^''z\'/''dz 
Jr 

Therefore, using this the former inequahty turns to be 



r — r 



/II 2 



d^r d^r' 



-s/2 



-k 

oo 

<MoY, 

-k 



,3n 



c- I / t^-'ds 1 c"(^-"-^) 

5c 



4-s 

<Mo(5,5,c)5]c"(^-^-^). 
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So the s-potential will be bounded whenever A — s — H > 0, that is, s < A — H and, 

4 -//< dim^^ graphs^ (1.110) 

Then comparing this equation against (1.109) it is: dim^j graph — A — H. Finally, 
from equation (1.107) we find that 

dimHe~\c) =3- H. (1.111) 

This expression exactly matches the one calculated by Constantin et al. (1991). Since, 
it replicates the equation (1.93) with (i = H as we proved above. 

We have shown here that the isotropic fBm not only provides stationary increments 
and reproduces the structure function for the stochastic refractive index but also gives 
the right fractal dimension associated to the passive scalars. Moreover, the structure 
function (1.99) gives a covariance function v which corresponds to a non-differentiable 
process. A simple calculation 

A^KIIr - r'll) = 2HA4v - r'|p^-2(a;, - x'^, 

shows this derivative is undefined whenever r = r', and because 

^l(llr-r'll) = ([/.(r + r') -/.(r')]^) =^(r,r)+^(r',r') -2^;(r,r'). 

The Cramer and Leadbetter's Lemma (1967) proves the refractive index is non-dif- 
ferentiable. This is the same with the isotropic fractional Brownian motion as it is 
proved in the Appendix A. 

Therefore, these reasons are enough to use this model as source in the Optics' 
differential equations. In particular, we will use plenty of it in the last chapter. 



Chapter 2 



Classical Methods Applied to 
Turbulent Lightwave Propagation 



The problem of light passing through a hollow made on a surface is explained in 
every Optics treatise. In this chapter we are going to make a brief introduction to it; 
afterwards, we will show how it is related to the Feynman's Path Integrals formalism. 
It is this technique which had proven fundamental studying image formation in the 
case of light propagating through a turbulent medium. Further ahead we will describe 
such a problem, and eventually use it to study some characteristic properties of the 
refractive index. 

2.1 From the Green's Theorem to the Feynman's 
Path Integral 

Assuming the polarization effects are negligible, the problem we have introduced is 
mathematically described as follows; 



AG + k'^G = 0, and 

G = all over the surface <t, 
G ^ M as r ^ 0, 
r (VG • cr - ikG) ^ as r ^ 0, 



(2.1) 
(2.2) 
(2.3) 
(2.4) 
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Figure 2.1: The surface where the boundary condition u is set is represented by its 
normal vector a. It is contained in the (x, 2/)-plane, while the 2;-axis is the direction of 
propagation. 



where G is the Green's solution to the wave equation (2.1), u is the boundary condition 
given by the hollowed surface (Figure 2.1). Therefore, the solution to the Kirchhoff- 
Huygens equation is written in this context as 

Anup = - d(T-{uVG), (2.5) 

J a 

where p is the point where we evaluate the propagated initial field. Now, let cr be a 
plane surface perpendicular to the 2;- axis, and by u{r, z) design the propagated field at 
a distance z from the initial (boundary condition) field u := u(r,0). So solution (2.5) 
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yields: 



1 f dG 
= ^ J(fpu{p,0) — 



z=0 



Then, we just need to find the Green's function. Applying the image principle to 
the free-space Green's function we build ours: 



G(p, z'; r, z) 



n r2 

where r\ — {p — r)^ + [z' — zY and r\ — {p — r)^ + {z' + z^. Thus, 



dG_ 

dz 



= 2ik 



ze 



ikr 



z=0 



rl2 



^ kr' 



where we have set r' — ri — r2 — y^{p — r)^ + z^. Now with this expression at hand, 
we may evaluate the solution along the 2;-axis, that is. 



ze 



ikr' 



r.12 



1 + 



2kr' 



(2.6) 



We can turn the pupil a into a function, and add it to the boundary condition u. 
Therefore, the integration is taken over the whole plane, which in turn can be expressed 
as the union of the sets {||p|| < z} and {||p|| > z}. Also, we are going to introduce 
two mayor assumptions: that u is symmetric around the 2;-axis, so da — 27rp dp] and 

kz ^ 1, i.e., the wavelength is smaller (~ 10~^ 10~^ m) than the distance to 

the origin. Under these conditions the integral over the first region has the following 
expression: 



Akr' 



kr' 



(2.7) 



Now, let us make the change of variables kp — x and set e = kz] also. 



kr' — e 



l + (^ 



1/2 



(2.8) 
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Thus equation (2.7) is rewritten as 



, , L (x\ expze 

H\\p\\<^}i^) = - y/^^^i^j — (I 



x\ expie{l + x^/e^^ 



1 + 



e(l + xVe2) 



(2.9) 



We may replace by a Taylor series in x/e all the expressions matching (2.8), and finally 
keep the terms up to the second order. That is, 



6*^ / X\ I X'^ 



/ x"^ 



— / xdx u(-) exp t { — = e / ( 



d^pu{p)expi i^P^ )■ 



Let us analyze the second case ||p|| > z: its contribution to the solution can be 
split. 



H\\p\\>^}i^) = -'^^^ I pdpu{p) 



Akr' 



r.12 



Akr' 



kr' 



— —ikz I pdpu{p) — — + kz I pdpu{p) 



Akr' 



Again, with the change kp = x both integrals are written as 

HM>^} - -^ri xdx u (- j j^^^ + r, xdx u (- j 



kr'^ 



,i(772+a;2) 



{rf + x"^) 2 



where 77 = kz. Instead, another change can be used here s — {rf + x^) 2 , and so 
ds — xdx/{rf' + x^)^. 



^{l|p||>4 



-17^ r u(k-\s'-r)'fA^ + v r u(k-\s'-r)'yA^. (2.10) 



This contribution is bounded; moreover, it tends to zero as 77 goes to infinity. With 
the condition A;^; ^ 1, we can make use of the mean value theorem and write. 



Us 



\p\\>z 



j ~ u{r]/k)7] < —i 



e^^ds 



+ 



V2v 



e^^ds 



V2v 



(2.11) 
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We observe that the integrals within the braces are compositions of Exponential- 
Integrals and linear functions, 



Ei{i\/2r]) and 
iEi{iV2rj); 



therefore, it is ii{||p||>z} ~ 0. 

Thus, the solution is just the contribution of the initial field enclosed by a sphere 
of radius z: 



u{z) 



k 



2Tiiz 



P\\<^ 



k 

d^pu{p, 0)exp ( i—p^ 



Furthermore, this integral can be extended to the whole plane assuming u{p, 0) ~ uop~^ 
as p — > oo. Remember, the solution we have build is valid whenever kz ':$> 1 is satisfied. 
Now, to evaluate the solution off the 2;-axis we just need to make the change. 



u{r, z) = 



k 



27:1 z 



d^p u{p,0) exp 



(2.12) 



This is the paraxial or Fresnel approximation, and it describes the free-space diffraction — 
from now on the term free-space will refer to a space free from inhomogencitics. Clearly, 
the phase term e*'^^ do not add information to the irradiance distribution — it is just the 
plane wave factor — so we will drop it off the paraxial approximation. Besides, setting 
u{r, 0) = e~**^^(5(2)(r) as initial condition let us build 



G{r, p; z) 



k 



27:1 z 



— exp 



.k 

'Yz^P 



that is, the Green function associated to the parabolic or diffusion equation: 



(^2tk^^+Vl^ G{v,p;z)^5^'\v-p) 



Also, the solution (2.12) has an alternative operator form, similar to the introduced in 
Section 1.2.2, 

ii(r,z) =exp(-i^V2)«(r,0). (2.13) 
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It can be deduced using the Fourier representation 



1 

2^ 



(2.14) 



Hence, 



u{r, z) 



471-2 

oo 

— Y- 



^—Y^ (-^y I d^K k'^M'^) e--"- 

^[d^.ui^)e.,(-i. 
271" Jr2 V 



r — t—K 

2k 



271 

I d^p u{p) [ d^n exp -iAv • (r - p) - i^fi^ 

47r^ ./ffi2 Jm L 2A; 



1 /" 2 f \ ^ 2nk 
a p u[p) 



47r2 



27riz 



iz 



I d^p u{p) exp 



exp 



like we said. 

This Green function will allows us to build a solution for 



d 

2ik— + Vl + kh{r,z) 
oz 



U{r,z) = 0. 



(2.15) 



This differential equation is the main subject of the remaining of this chapter. It is 
obtained from changing the vacuum refractive index in equation (2.1) by an inhomoge- 
neous one, like those we described in the last chapter; afterwards, the laplacian operator 
is approximated under the condition kz ^ 1 and the multiplicative term results from 
writing the solution as e^^^U . 

Now, we will follow Charnotski et al. (1993) building the solution to (2.15) from 
the free-space Green function. Suppose we want to find the scalar wave field at a point 
(R, L) given the boundary (initial) condition m(p, 0) = uo{p)- We will construct the 
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field at that point by subdividing the interval [0, L] in subintervals [zj, ^j+i) of length 
= L/N, with N a large integer. Let us use equation (2.15) to estimate how the 
field propagates from a point Zj to the next -Zj+i, 

u{r, = u{r, Zj) - ^ [V^ + ^'^(r, z^)] u{v, z^) + ^(Az^) 



iAz 



exp i ^ + A;^e(r, ^j)] |> ^(r, Zj); 



(2.16) 



this last expression is exact. Remembering that e^^^ = e^^'^^e^e^ we can detach 
both operators above; moreover, because both operators in above equation are of order 
one in Az the commutator is of 0{Az^), and so its contribution can be neglected. 
Therefore, applying the homogeneous paraxial solution (2.12) and property (2.13) we 
turn (2.16) into: 



u(r, Zj+i) = exp 
= exp 



/ 



' kAz 

■ kAz 

I 

2 

k(fp 



e(r, Zj) 
e(r, Zj) 



exp(^z^V^^ u{v,Zj) + 0{Az^) 



2m z 



M2 



(fp u{p, Zj) exp 



k 

i — (r - p)^ 
2z^ 



, , ik 
. Mip, 2;, )exp — 
27riA^ ^2 



Az 



+ Az e(r, Zj) 



+ 0{Az'') 
+ 0{Az^). (2.17) 



We have build a recursive algorithm to estimate the field at the arriving point Zj from 
its predecessor at Zj+i. The field at the point (R, I>) should come after N iterations 
from the initial field at 2:0 = 0, that is. 



m(R,L) 



k dPpN- 



k <PpN-2 



k (Ppo 



]K2 27iiAz 2niAz 2niAz 

X exp i 77^7 [(^ - PN-if + (Piv-i - Piv-2)' + • • • + (Pi - Po)'] 



2Az 



ikAz 



+ 



+ 



[e(R, zn-i) + e(pjv-i, ^iv-2) H h e(pi, 0)] \ lio(Po) 



(2.18) 



NAz)' 
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Afterwards, for any function given p{z) the set {p(0) = Pq, . . . , p{zj) = pj, . . . , p{zn 
L) = R} will define an interpolating function (Figure 2.2), with constant derivatives 
within the subintervals [zj,Zj+i), which converges to it as ^ oo. Also, note that 
with A^ growing we have ^/NAz 0, and then the terms inside the argument of the 
exponential in the latter equation behave: 



lim -— 7 (p,+i 

j=0 



lim > 

N-^oo 



j=0 



Az 



Az 



dz 



dp{z) 

dz 



n 2 



AT 



lim Az > eip. 



dz 



e{p{z),z) 



(2.19) 
(2.20) 



Nevertheless, these convergences are not enough to prove the existence of a limit for 
the expression in (2.18) when A^ — ^ oo. We may think that the following definition 
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plays the role of a measure 



N-l , ,2 

V'p{z) := lim n (2.21) 

k=0 



but it is not true. In fact, it is exp |(iA;/2) J^dz [dp{z)/dz]^^ V'^p{z): 
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or, 27riL 
°° kN d^uj 



oo ./ — oo 



27riL 
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ml 



j=0 

N-l 



kN Ujduj 



iL 



lim Uil)^!, 

AT— »oo 



(2.22) 



j=Q 



here we have made the change p^ — H — Yl^=i — the notation for the integral on the 
left-hand side stresses the fact that the starting point is not fixed as one can see from 
the definition (2.21). 

We conclude that the hmit to the expression at the right side of (2.18) exists, and 
we write it as: 



u 



p{L)=-R 

v'p(z) uo[pm 



ik , \dpiz) 
X exp< — I dz — i-^ 
T ./n dz 



+ y ^ dze{p{z),z) 



(2.23) 



This is the solution to the wave equation with a non-constant refractive index in the 
paraxial approximation, (2.15). Even we can retrieve its associated Green function by 
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making t^o[p(0)] = S^^^ [p(0) - Rq] , i.e., 

-p(L)=R 



G'(R,L;Ro,0) 



ik 

X exp — J dz 



n2 



dpjz) 
dz 



y io e(p(^)'^) 



Now, the reciprocity theorem (Sommerfeld, 1949) states: G(R, L; Rq, 0) = G(Ro, 0; R, L). 
Applying it to the former equation yields, 

G(Ro, 0; R, L) = / V^p{z) 6^^^ [p{L) - R] 

, n 2 



, ik rcip(;^) 
X exp < — i dz — i-^ 
^ 2 7o [ dz 



+ y ^ dze{p{z),z) } . 



We observe here that the 5-function alternatively replaces the conditions of the form 
p{z) = R, so symmetry imposes 

^(R, L; Ro, 0)=j Vpiz) 5^^\p{Q) - Rq] 5^^\p{L) - R] 





'dp{z)' 




dz 



+ 



ik r 



dz e{p{z), z) 



In this context, the normalization can be rewritten in one of the following forms: 

fp(L)=R 





'dp{z)' 




dz 



-r 

= J V'p{z) [p{L) - R] exp I ^ J^z 

= 1 1)V(^) 5^'\p{0) - Ro] exp I ^ ^^rf;. 



dp{z) 
dz 



(2.24) 
(2.25) 
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Afterwards, the general solution to (2.15) is 



Next, we are going to show how this construction help us determine the irradiance 
pattern over a screen. 



2.2 Image Formation using the Feynman's Path In- 
tegral Representation 



In the following we are going to change the representation space we are actually using 
by another called velocity representation space. This representation was introduced 
and used frequently by Russian scientists — it is just a functional change of variables. 
It has proven extremely useful in handling propagation problems (Klyatskin, 1970, 
1975). Now, let us introduce it: 




(2.26) 



' dp 

dz 




R, p(0) = Ro, and 



(2.27) 



I v(0) 



Also, this new variable requires its own normalization, that is. 




with V^viz) = \miN-.ooY{^=i{kLd'^Vj/2mN). 
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This change transforms the inhomogeneous Green function (2.24) into; 
G(R,L;Ro,0) =exp|^(R-Ro)='| J V\{z) 5^^^ Jo"^^^^^^ 

J z 



X exp <^—fdz v^{z) + — f dz € 
2 Jo 2 7o 



(2.28) 



This representation allows us distinguish between the contribution made by the inho- 
mogeneities from the free-space propagating wave. Now, let us calculate the free-space 
Green function to introduce the procedure we will afterwards follow in more complex 
situations. As always, we start dividing the interval [0, L] in N new subintervals, but 
this time we will also introduce the following Fourier transform. 



2n 5(2) 
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We have then. 
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So it is appropriate to write the Green function as: 



^(R, L; Ro, 0) = ^(Ro, R, L) G'o(R, L; Rq, 0), 



(2.30) 



where Gq is the free-space Green function and 



^(Ro,R,L) -.^^-^ j V\{z) 5^'^ 



I dz y{z) 
Jo 



X 



ik 



ik 



X exp < — / dz V (z) -\ / dz e 



z (L — z) 

-R + ^ ^ ^ Ro + / dr]v{r]),z 



; (2.31) 



therefore, we have concentrated aU the random features of the medium within this last 
function — it is set equal to 1 when perturbations are absent. 

Before going further, we have to describe the environment in Optics where ap- 
proximation (2.15) holds. As in Section 1.3.1, we write the permitivity as a constant 
term — we will assume equal to one — plus another term e containing all the information 
coming from the medium. There, we consider the propagation of quasi-monochromatic 
lightwave radiation with frequency cu, that is, 

{A+[l + e{r,z)]k^}S{r,z) ^0, 

where S{r, z) is the scalar electromagnetic field — there is no polarization here. Suppos- 
ing that backscattering is negligible, and thus the propagation has a preferred direction, 
let us say £(r, z) = E{r, z) e*('=^~'^*), the latter equation changes to 



[A + A;=^e(r, ^)] E(r,^) =0. 



(2.32) 



Hence, under the condition kz ^ 1 we retrieve equation (2.15). Thereafter this is 
a stochastic parabolic equation, and it can be solved when e is a markovian^ process 
along the propagation axis as it was shown by Rytov et al. (1989). Moreover, its 
solution coincides with the deterministic solution we have just shown. Charnotski 
(1991) extensively discusses the apphcations of the path integral formulation to this 
problem. In the Appendix B we summarize how he characterizes the cohabitation of two 
regimes: weak and strong. Basically the technical differences between both regimes are 
^See Chapter 3 for a detailed description. 
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the following. The weak regime cumulates a variety of methods grouped under the tag 
of Rytov's formalism — whose central idea is to solve equation (2.32), or an equivalent 
form, using a Taylor-like series expansion — in practical terms this distinction provides 
a way to select the best method for solving a particular problem. While the strong 
regime has, until now, one tool: the path-integral approach. Nevertheless, it is this 
approach the only one covering both regimes. 

Finally, let us build the irradiance distribution from the solution to the inhomo- 
geneous wave equation (2.32). According to definition (2.30) and equation (2.26) we 
have: 

E(R, L) = / (fRo E{Ro, 0) ^(R, L; Rq, 0) 

= / (fRo E{Ro, 0) Go(R, L; Rq, 0) ^(Rq, R, L), 
and so the intensity function is 

/(R,L) = E*{R,L)E{R,L) 

^ [ [ (fR'o d^Ro E*(R'o, Q)E(Ro, 0) ^^(R, L; 0) ^(R, L; Rq, 0). (2.33) 

Now suppose the coherence time is much smaller than the characteristic time scale 
T of the detector, i.e., u!~^ <C T. Thus the irradiance pattern observed is time averaged; 
furthermore, it can be assumed ergodic. So, if the source is spatially incoherent, the 
mutual intensity is written (Goodman, 1985) 

^*(R'o,0)E(Ro,0) = A(Ro) S^^\K - R-o), (2.34) 

the overbar means the ensemble average. Besides, the function A{r) has dimensions 
irradiance Xcirea, that is, is an intensity distribution; moreover, in most of the cases 
it is proportional to the initial irradiance distribution at the pupil: 

A{R) = -/o(R). (2.35) 

TT 

Nevertheless, for the current discussion we will keep the distribution function A, and 
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so we use the property (2.34) to write the intensity function as 



with G as in (2.30). 

In the next section we will apply these results to the case of image formation in 
self-image system. 



2.3 Intensity Distribution of Self-image Systems into 
Turbulent Media 



Here we will describe light propagation through a Lau-like arrangement, that is, two 
Ronchi grids out of phase half a period within a turbulent medium. Also, we will 
inspect how degradation produced by the turbulence can be estimated in terms of the 
spacing between two parallel grids, the number of lines per millimeter, and C^, the 
structure constant of the medium (Perez and Garavaglia, 1999). 

Now, we introduce the optical system for the present discussion: as it is sketched 
in Figure 2.3, it consists of a Lau system of grids — we have no lens here — separated by 
a distance L and half a period out of phase. That is, each grid can be thought as the 
negative image of the other. Finally, between the grids there is a turbulent medium 
with structure constant C|, roughly homogeneous in the plane perpendicular to the z 
direction. Ai z = I there is a screen where we want to observe the system behavior. 
We suppose that the medium between the second grid and the screen plane is free from 
turbulence. 

In this problem there are only two optical elements: the grids. Both of them have 
the same tramittance function, it can be modeled as a family of square wave functions 



the parameter a describes the relative phase and 2d is the grid period. There is a 
difference between both grids when we turn them into mathematical objects for our 
problem: the light passing through the first grid coming from a spatial incoherent source 




(2.36) 




(2.37) 
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Figure 2.3: The optical system employed in this work; two grids separated a distance L 
with equal amplitude tramittance functions {Ai is the outgoing intensity distribution 
and A2 the tramittance function, these both have a period 2d) and / is the intensity 
distribution at a screen located at z — I. 

can be expressed as a delta correlated coherent function with intensity distribution 
Ai{v) given by 

A(r) = ^oto(e. • r)5^xD(r), (2.38) 



where 

I U otherwise 

Aq is the maximum value for the irradiance distribution and D is the size of the 
rectangular grid; while the tramittance A2 of the second grid is modeled using the 
same functional relationship, but with Aq — 1 and setting the phase a — d. 

The irradiance function /(R), as we showed arriving to equation (2.36), can be 
expressed as follows: 

7(R)= / d\Ai{r)\Gtot(R,l;r,-L)\''. (2.40) 



This time the total Green function Gtoti^, ^; r, — L) comes from combining the free- 
space Green function G'o(R, I; r', 0), which corresponds to the zone between the second 
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grid and the screen, and the turbulent Green function ^(r', 0;R, —L): 

Gtot(R, I; r, -L) = / GoiR, I; r', 0) A2{r') G{r', 0; r, -L). (2.41) 

Remember that the irradiance distribution Ai and the tramittance function A2 were 
built from delta correlated fields as was discussed in the last section. Now, combining 
(2.40), (2.41) and the turbulent Green function definition (2.30) yields: 



{ZTT) ( 1j JdxD Jm?JDxD 



(27r)2/2L2 



• I R' + - 



exp 



exp 



(2.42) 



x^(R'--,r,Lj(7*(R' + -,r,L), 



here / stands for the relation 1// = l/l + 1/L. This equation is similar to that found 
by Charnotski (1996) but here we have not lost the phase term exp(— i^r' • R'), just 
because we are working with grids instead of lenses. 



2.3.1 The non- turbulent case 



Before treating the turbulent problem we are going to consider light propagation in the 
absence of turbulence. Our goal here is to inspect the role of each grid in the image 
formation process. Let us assume that 5'(r, r', \z — z'\) = 1, thus equation (2.42) takes 
the form: 



I^ztt; ( L, Jdxd Jdxd 



, (2.43) 



where 



Ct(r')= / d^R'td 
Jdxd 



e, ■ ( R' - ^ 



■ I R' + 2 



exp 



i|r'-R^ (2.44) 



is a complex correlation function. 
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Now, let us carefully inspect the former equation. Noticing that 



/ ' 

J DxD 



(frto{ex •r)exp 



—tr ■ — 



kr' 



■r') ^ 2.^4(4 r'), (2.45) 



where ©(^^(r) = 0{rx)'d{ry) is the two dimensional Heaviside function and D = -De^ + 
Dey. We shall rewrite equation (2.43) as, 



27r/fR) = 

= / 



An k^ 



therefore, with the change of variables K = kr' /I the function between brackets results 
to be the Fourier transform of the irradiance /. That is. 

Let us inspect each term in this transform. From the definition (2.37) we observe that 
the Fourier transform for the transmission function to rests on the decomposition 

DxD 
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(2n+ l)7r 
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-X 



2iexp 



.(2n + l)7r 



-X 



+ exp 



(2n + l)7r 

-I ; X 

a 



(2.47) 



As far as we are concerned with the effect produced by the grids, we are going to 
neglect the low spatial frequencies related to finite size effects, that is, we will take the 
limit 1/\k\D — > each time we calculate the resultant irradiance — equivalent to the 
condition D ^ X. Otherwise, the optical system acts as a filter transmitting only a 
discrete numerable set of frequencies. 
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Afterwards, using the dimensionless variable r]/X — Kwe finally have, 









exp 







n=—oo 

The last series comes from the limit D/A ^ oo, its terms are: 



(2.48) 
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(2.49) 
(2.50) 



Similarly, the complex correlation can be evaluated. First, we multiplicate the 
tramittance functions td to obtain: 
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= Mo+ 5^ Mf)+ ^2"'"^- (2-51) 



rn,n=— 00 



These terms are directly Fourier-transformed because of the definition (2.44) under the 
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condition we have given above, i.e.: 
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(2.53) 
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(2.54) 



Finally, the Fourier transform of the irradiance distribution (2.46) is obtained multi- 
plying the tramittance and complex correlation transforms. Nevertheless, we should be 
cautious. The approximation we suggested induce the constant terms to produce more 
delta functions than the Fourier integral is able to handle. Therefore, we will redefine 
those problematic terms. Let us start with y-axis, because all tramittance functions 
depends on the x-axis we can write. 



Vly{Ry) = 



IL 
IL 



[ 2 ' 2 



^drydR'ydr'y exp 



{A 
— / dCexp 



.RyTy \ > 

I ^ L f 



drydR[, 



k^/2'KAo 
E 

kfV2^o 
IL 



RL 



/;;../<.(^.^^)e(f-«;)e(f.^. 



We observe from above that eliminating the effects from the edges corresponds to the 
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condition —D/2 < Ry < D/2, there we have the constant intensity distribution: 

kfDy/2¥A^ 



IL 



Furthermore, a constant term will appear when multiplying both first terms in the 
tramittance and correlation functions, (2.49) and (2.52), 

— this is the contribution from the edges to the irradiance. 

Also, it is worth noting that the products if^ Mi and M^i^^^ are zero because (2n + 
1) 7^ for all n e Z, but the remaining cross product by constants contribute to the 
irradiance distribution on the x-axis: 
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(2.55) 



We changed the sum in m to the integral, inw — (27rA//d Z)m, because of the condition 
d/X^ 1, where d is defined through N{2d) — D — N is the number of hues of the grid. 
We find two others non-zero contributions to the irradiance, following the same 
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procedure as above: 
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(2.56) 



The irradiance is real function, so the exponential in the latter equation should be real; 
it only happens when 
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(2.57) 



for s G Z"*", any other choice in the quotient will make the series vanish. Finally, we 
write equation (2.56) as 
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. (2.58) 



2.3 Intensity Distribution of Self-image Systems into Turbulent Media 65 



And the second non-zero term is, 



(2) 



nl 



J2 4"^Mi"'"'^ 



n,m,m'=—oo 

oo 



E 



nl ^ (2n+l)(2m+ll 

n,m=— oo ^ ' ^ ' 



Vx + 



{2n + l)7iXL 
dl 



X 



E 



m'=— oo 



(2m' + 1) 



exp 



.2(m — m'){m + m! + l)7TXf 



X S 



(2n + l)7rAL 2(m + m' + l)7rA/ 



dl 



dVh 
nl 



E 



n,m=— oo 



-1 



\m+n+l 



(2n+ l)(2m+ 1) 



dl 



Vx + 



X 



/_ 



dw 



cxp < t 



X exp 

dVh ^ 



/2 



2nXf 



2nXf 



dl 



m — w 



(2n + l)7rAL 
dl 

/27rA/\ 



\ dl ) 



{m + 1) + w 



nl 



E 



d<^w — 

j^\m+n+l 



(2n+l)7rAL 2(m+l)7rA/ 
dl dl 
(2n + l)7rAL 



n,m=— oo 



1 



X 



(2n + l)^-(2m + l) 



(2n+ l)(2m+ 1) 
expi 27r(2A; + 1) 



Vx + 



dl 



xexp<. 



(2^ + l)2^-(^ + l) 



/vrAL 



(2m + 1) - (2n + 1) 



L 
2/ 



(2.59) 



and these exponentials should give us a real number. It takes some algebra to rewrite 

them as _ ^ 

" AL _ \ f L\ XL 

1?' ) \2f) + d^. 



exp < ^TT 



If we now assume that (2.57) is fulfilled then it happens that only XL/d"^ — 2p + 1 
makes this exponential real. On the other hand, it also makes the series in m from the 
former equation to have a singular term. Nevertheless, we realize that the condition 



(2m+l)(2s + l) = (2n + l) 
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— according to (2.57) it is L// = (2s + 1) — gives a zero term in the original series. 
Thus, we finally have the series 
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This series can be reduced: because 
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and knowing (Gradshteyn and Ryzhik, 1995) that 
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Therefore, the equation (2.59) is 



(2n+l)= 



2s(2n+ l)7rA 
d 



(2.61) 



under the conditions L — (2p + l)ri^/A and (2s)/ = L. 
Now we can change the relation (2.57) to 



I = 



(2s + 1)' 



(2.62) 



which automatically makes the equation (2.58) vanish. That does not happen with 
(2.59); furthermore, it induces a relation of the type 

L = —2p (2.63) 
A 

with p e Z. It also forces the distance between the last grid and the screen to be 
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(2.64) 



Observe that q and p must be simultaneously odd or even. Then the series (2.59) is 
written. 
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Studying again the sum in m: 
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(2.65) 



Finally, we are ready to write the complete Fourier transform of the irradiance 
distribution due to the grids Igi 

Case L = (2s + 1)/: All terms are zero but the constant (2.55). Supposing the rela- 
tions (2.63) and (2.64), imposed by equation (2.59) still applies, then 
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Case L — {2s)l: The full irradiance Fourier transform has two possible expressions. 
When L = (2p + 1)£/X, 
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(2.67) 



The other non- vanishing contribution has a simpler expression. 
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for any other L. Moreover, the above expression can be simphfied by tuning L, we 
choose it in a similar fashion as in equation (2.63): 



^ 2d\ ^ , 2d\ 
L — — r-2p and / = ~t~^Q.i 
A A 



(2.69) 



which maximizes all the series terms. 

Now, we can recover the full irradiance distribution and calculate the visibility in 
each one of the exposed cases. Most of the situations will give us a constant intensity 
distribution, that is. 



7(R) 



ae + /.)(R) 

7e(R) = 2 



"(2s+ l)7r" 




. (s + l)p . 





2N + 



>+i: 

80F 



{2s + l)7r 
(s + l)p 



iV/o, for L = (2s + 1)1 
for L ^ (2s) L 



Thus, we obtain a non-constant irradiance distribution only with the condition L — 
{2s)l. We have seen there are two possible solutions: if L = (2p + l)d^/A is 
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, (2.70) 



or expressed in terms simple periodic functions as 
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for -d/4s < Rx < d/As, and 
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for d/As < Rx < d/2s or —d/2s < Rx < —d/As, and so on. 

On the other hand, when L — (2s) Z we choose the relations (2.69) and thus 
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(2.71) 



Again, we express it using periodic functions 
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for —d/2s < Rx < d/2s. This irradiance pattern is a IG*'^ part of the latter (Figure 
2.4). 

The quality of an irradiance pattern — the contrast — produced by a system of grids 
is quantitatively measured the visibility V defined by Michelson (Hecht and Zajac, 
1986): 
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(2.72) 



where /max and /min are two consecutive maximum and minimum. It is equal to zero 
for all distances but those described above. When L is an even number of d^/X we 
have 
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(2.73) 



the with of the grid d, the wavelength A and the distance to the screen I completely 
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t 




shape is found whenever L = {2p)2(P/\ and / = {2q)2d'^/X, but the shape has less 
contrast — it is 1/16 of the latter, c) The triangular shaped and parabolic teethed 
functions contribute to a) but just the former to b). The parabolic teethed function 
weights considerably less than the triangular so its contribution is almost negligible. 
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define s, and thus the visibihty. Otherwise, when L = {2p)2(P/X the visibihty is just 

V = s. (2.74) 

Therefore, it only depends on the quotient between L and I. The nearer the screen to 
the last grid the higher is the value of the visibility 

2.3.2 The turbulent case 

The statistical averages we will use here are understood as long exposure time-averages 
(Roddier, 1981). Furthermore, the only relevant assumptions about the permitivity e 
is being Gaussian process and markovian on the 2;-axis (Tatarski and Zavorotny, 1980). 
From equation (2.42) and the definition (2.31) we must evaluate the following^ 
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(2.75) 



Now, the averaged terms within the exponential can be rewritten using the markovian 
property, that is. 
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^For any given Gaussian process X, we have (expiX) = exp— (X^)/2. 
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and 
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where A{p, z) is defined as in Appendix A, but with a 2;-axis dependence. Thus, we 
write the exponential term as 
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Therefore, we can introduce the linear change of variables: Vi — V2 = v and Vi + V2 = 
-V. Because t^K^) - vl{z) = (v2 - Vi) • (vi + V2){z) = 2v{z) ■ y{z), we turn (2.75) 
into 
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We can group all the dependencies on V and integrate. As in the classical calculus it 
give us a delta function, that is. 
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the delta for functionals. Moreover, when one of the extremes is fixed, as in our case, 
the path-integration of it is not equal to one: 
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henceforth, A^o = N^/^-n Lr and so it is 
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The latter property and the delta's definition allow equation (2.75) achieves its final 
form, 



(^g(R' - ^, r, g*(R' + ^, r, ^ = exp / ^ [(^ " ' 

= exp-D(r',L)/2 (2.79) 

We shall proceed to evaluate the mean irradiancc function. We have shown the 
average adds a function dependant on the coordinate r', then we arrive to an equation 
similar to (2.43) but with an extra term: 
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Its Fourier transform is now straightforward, 
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Whether it is the case of equation (2.67) or (2.68) the exponential contribute to each 
term of them with 

exp -D{Q{2n + l)de^, L) /2, 

here Q is an integer satisfying one of the conditions we have given. Assuming the 
structure constant Cl{z) is roughly homogeneous we can write 

D{re,,L) ^ ^^(^^^) sm^ k^C^Lr''+' = D.r'^+K (2.82) 
^ ' ^ 4r[(// + 3)/2]' 2 ' ^ ^ 
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Henceforth, for L = {2s)l is 
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with Q = {2p + 1) when L = {2p + l)d'^/X, and when the relations (2.69) are satisfied 
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(2.84) 



with Q — Ap. 

Thus, only a finite number of terms contribute significantly to the image forma- 
tion. The exponential term in both series, (2.83) and (2.84), plays the role of a cutoff 
smoothing the original irradiance pattern. The integer N ~ {10/0^)7^ {2Qd)''^ is a 
good measure of this cutoff — terms beyond that number add corrections of order less 
than 10~^ to the actual value. Also, it makes the irradiance extremely sensitive to 
changes in the wavelength and the structure constant. Figure 2.5 shows the difference 
between the patterns generated by infrared and ultraviolet wavelengths for the same 
geometric arrangement. 

Afterwards, we can estimate the visibility. The visibility in the turbulent case is 
smaller than in the non-turbulent one because of the cutoff, and it turns smaller as the 
wavelength decreases: for L— {2p+ l)d'^/X, 



V = 0.69742, Vi.2;.m = 0.67865, and V4oonm = 0.63036; 
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Figure 2.5: The figure compares the irradiance patterns for two very different wave- 
lengths, 400nm (soft-ultraviolet) and 1.2/im (red), given a fixed geometric configura- 
tion: L = 0.976m and d = 0.625 x 10~^m. The red-wavelength (p = 1 to reach the 
distance L) function has been mirrored to compare against the other two. 

and for L = {2p)2<f/\ 

V = 1, Vi.2;.m = 0.88163 and V4oonm = 0.66679. 

Amazingly, it is the second irradiance distribution pattern (2.84), which has a flattened 
pattern, more sensitive to changes in the wavelength and turbulence behavior against 
what their patterns suggest. 

The cutoff also depends on the geometry of the system. Two instances are relevant; 
as — i> oo the visibility goes to zero, otherwise if c/ — > it takes the same value as 
in the non-turbulent case — equations (2.73) and (2.74). These results show us how 
the geometry influences image formation in a turbulent media. The behavior of the 
visibility is in agreement with the results of Zavorotny (1988) for an infinitely extended 
source as it vanishes when d goes to infinity. Moreover, if d is small enough the effects 
of the turbulent medium vanish and the visibility recovers the value it takes in the 
absence of turbulence. 

Finally, here we have established the conditions for image formation in a Lau-like 
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Figure 2.6: The first graphic displays the irradiance patterns for a A = 400nm wave- 
length and the second for A = 1.2/im. The degradation is clearly observed in the first 
example, but hardly can be seen in the red wave length. 
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arrangement. For a visibility different from zero, the separation between grids, L, and 
the distance from the last of them to the screen, Z, are related by the condition (2.57). 
We observe the appearance of a characteristic length ^ j A, it is called Talbot distance 
and is widely present in grids systems. Only on integer multipliers of it we have found 
a non-zero visibility. In these situations we were able to express the degradation in 
terms of a few variables: the physical and A, and the geometrical L and d. 

Also, the mean irradiance is exact: either it is useful in both strong and weak 
regimes. Equations hke (2.83) and (2.84) provide us with a new way to calculate the 
structure constant of the medium at laboratory from a density section of an image. 
Indoor experiments carried out with laser beams through turbulent medium (Consortini 
et al., 1990, 1996) are based in measures of their wander and thus an statistical analysis. 
While ours just needs an interpolating Fourier polynomial. 

We have given an introduction to the classical methods in turbulent propagation 
based on a markovian model. In the forthcoming chapters we will introduce processes 
with memories to accurately resemble the model we introduced in the first chapter. 



Chapter 3 
Stochastic Calculus 



We have shown that defined the turbulent refractive index as a member of the family 
of fractional Brownian motions it is not differentiable. Furthermore, we usually find in 
Optics derivatives of the refractive index within differential equations, but when the 
media is turbulent these equations are undefined in terms of the Classical Calculus. 

For instance, let us suppose it is possible to define the derivative of a fractional or 
standard Brownian motion, the noise: . Thus, the integral equation associated to 
x{t) = B^{t) x{t) is just 

x{t)=xo+ [ x{s) B"{s)ds = xo + [ x{s)dB^{s), 
Jo Jo 

for the last term above we have assumed that the change of variable formula is still 
valid. So, it is the existence of this kind of integrals what we should try to verify. If 
we attempt to define this integral as the limit of the Riemann series, 

n-l 

J2<m{B"{t^,)-B"{t^)), (3.1) 

i=0 

its existence can not be proven in general. 

Nevertheless, conditions over the argument function x(t) for the existence of this 
type of integrals are now well established, and a Stochastic Calculus can be build from 
it. This calculus and how it can be used to solve stochastic differential equations will 
be described next. 
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3.1 Introduction: White Noise and Brownian Mo- 



In this section we will present the stochastic analysis for the standard Brownian motion, 
and set the notions that will be later extended to the more general fractional Brownian 
case. 

The theory of distributions had provided us with derivatives for functions without 
them in the classical sense. Therefore, it is natural to propose the white noise as a 
distribution, but to do so we must also give the right abstract probability space. It 
was Hida (1980) who first used this idea as the building block for a stochastic analysis. 
Here we are going to build such a space and show how it allows define integrals in the 
sense of (3.1). 

Let 5(M°') be the Schwartz space of rapidly decreasing smooth (C°°(M'^)) real 
valued functions on and let us choose its dual <S*(]R'^) — the space of tempered 
distributions — as the probability space Q. We represent with (o;, 0) = uj{(j)) the action 
of the elements of the dual, u e »S*(M'^), on the functions belonging to 5(M'^). 

Of course, to properly define the probability space we have to attach a cr-algebra 
and a probability measure. The former is straightforward, we just use the family of 
Borel subsets B{S{W^)), and associated to this algebra we need to prove the existence of 
a measure. The Bochner-Minlos theorem (for a proof see Holden et al., 1996, Appendix 
A) shows that such a measure, exists; moreover, it has the following property: for 



where || ■ 1| is the norm in L^(R''). Therefore, we call the triplet {Q, B{Q), /j.) the 
1- dimensional white noise probability space. 

The probability measure is a Gaussian measure on iS(]R'^): we just have to evaluate 
the finite dimensional measures. So, let us take a set of functions .^i, ■ ■ ■ ,C,n G 5(M'^) 
such that they arc orthonormal in L^iW^). Now, given a random variable uj, we can 
project it into the finite random variable ((a;, ^i), (a;, ^2), • • • , (i-^, Cn))- For any smooth 



tion 



all (f) G 5(M''), 




(3.2) 
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function / e C°°(M'*) we have, 



E[/((-,6),---,(-,U)] = 




— we used above the properties of the Fourier transforms. Thus, we have found the 
n-dimensional Gaussian measure 

dA"(x) = (27r)""/2e-^"""'(ixi • • • dxn- (3.4) 

With the same procedure we can prove that if (p G L'^{M.'^) for any succession 
0„ G 5(]R'^) such that 0„ in L^(]R'^), then 3 hm„_>oo(^, 0n) '■= (^^,0) in L'^il^)- 

Let us introduce now the 1- dimensional (d-parameter) smoothed white noise. It is 
a map w : S{R'^) x S*{M.'^) M given by 

w{(f)) = w{(j),u;) = {uj,(j)); u E S*{R'^),(f) E S{R'^). (3.5) 

Now, we define the following process 

B(x) -.^ B(xi, . . . ,Xd,u;) ^ {u;,X[o,xi]x-x[o,xd]), (3-6) 

for X = {xi, . . . ,Xd) G M.'^, where x is the index function: gives 1 when x is inside 
the box [0, xi] X • • • X [0, Xd] and zero otherwise — when < it is convention to 
assume [0, Xi] represents [xi, 0]. This process has a continous version which turns to be 
a d-parametcr Brownian motion. 

It is evident from definition (3.6) that this process is almost surely zero at x = 0. 
Also, the process satisfy definitions (A. 5) and (A. 7) in their d-dimensional equivalent 
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form, that is, 



E 



5(x) 



and E 



a 

B(x)5(y)] = \[m.iTi{xi,yi}. (3.7) 



Checking these properties is straightforward, we choose x*^^^ , • • • , x*^") e M*^, and con- 
stants ci, . . . , c„ G M, so we build the index functions: x*^*^ — 'X^^Qxi]x---x{OxdV therefore, 
we compute the n-dimensional characteristic function. 



E< exp 



1=1 



E<^ exp 



i=i 

exp(^-iiiX:Qx«ir) 



exp 



n 



CiCj 



1 



exp ( "2^^^^ ) ' 



where c = (ci, . . . , c^) and V is the symmetric nonnegative definite matrix defined by 



Therefore, B is Gaussian with mean zero and covariance matrix given by V. It is 
better now, instead of directly evaluate the covariance, calculate the variance of its 
increments. So, making use of (3.3): 



E 



(^(x) - B{y)f 



= E[(-,X[0,x] -X[0,y])^] 

= ||X[o,x] -X[o,y]||^ / u'^d\^{u) 

d 

= ||X[0,x] -X[0,y]|P = H'^^ 



X[0,x] -X[0,y] \2 
-X[0,y]||^ . 



1=1 
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where [0, x] = [0, Xi] x • • • x [0, x^]. Thus, the Brownian motion covariance follows from 
the variance we have found. Finally, the continous extension to the process comes 
from the application of the well-known Kolmogorov's continuity theorem, and makes 
the continous version i?(x) the desired d-parameter Brownian motion. 

With this definition of Brownian motion we can define the Wiener-Ito integrals. We 
will simplify the following exposition setting d — 1. Let (f) e L^(]R) be deterministic 
with finite support set, let us say the interval [a, b]. Now, we build the succession: 

n 

0nW = $]0(t.)XM.+i)W' (3-8) 

i=l 

where a < ti < ■ ■ ■ < t„ < tn+i = 6 is a partition such that max |tj+i — tj| ^ as 
n oo. This family of functions belongs to L^(]R) and converges to there. The 
requirement for {0„}„eN being in iS(]R) is found making the edges of the function 
smooth in a neighborhood of the interval and approach to a step function as n grows. 
Let us omit that step to simplify the exposition, therefore, 

n 

{ujAu) = E'^i^o mi^i) - B{ti)) ^ 

i=l 

in — in mean square^. Thus, we can put 

(j){t)dB{t,u}) -.^ {u},(l)); ueS*{R), (3.9) 

The same arguments can be used with the d-parameter Brownian motion to define 
the stochastic integral in the same way. Moreover, we can integrate by parts — provided 
the pathwise integral coincides with the L^-stochastic integral — and get 

w{4>)^ [ (j>{x)dB{x,uj)^{-lf I {x)B{x,uj)d''x 
jRd J^d dxi--- dxd 

dxi---dxdj \dxi---dxd /' 

See footnote 2 at Chapter 1. 
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•, •) is the inner product in L^(M'^); thus, in the sense of distributions we write 



w 



dxi ■ ■ ■ dxd 



(3.10) 



Now, we would like to replace the deterministic function be a stochastic process 
f{u!,t). For most applications is enough to prove this replacement is possible for a 
closed set, say T — [0, 1], and we will do so. The extension, known as Ito integral, is 
possible whenever the process has the following properties: 

i) Given the set J-'t = {B{s) : < s < t}, then / is .7-^- measurable for any t 
{T^:= {f),0}). 

a) the map {io,t) — > f{uj,t) is B(R) x .Fr-mensurable. 

in) E[Jj,p{u},t)dt] < oo. 

This conditions are enough to guarantee the existence of the following limit, ltd integral, 

„ 2"-l 

/ f{u;,t)dB{t):= hm J] /(u;, j2-") [i?((j + 1)2"") - i?(j2-")] in L^(/x). (3.11) 

Irr n— >oo ' 



The choice of the step function ^^=7 /('^)i2~")x(j2-",(j+i)2-"]i which also converges 
in L^(f2 X M) to /, is crucial here. For it not only assures the limit in mean square but 
provides the isometry property, 



E 



f{uj,t) dB{u,t) 



E 



f{u,tfdt 



and also 



E 



f{u,t)dB{t) 



= 0. 



(3.12) 



(3.13) 



Alternatively it can be proved (Nualart, 1995) that the step function process 

r-i2-" 



2"-l 



" ' ' f{uj,s)ds ) Xo-2-",o-+i)2-n] 

0-1)2-" 



is also ^j2-"-adapted, converges to / and gives the same limit integral (3.11) with the 
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properties above. Other approximations to the process can be built, but they do not 
obey the latter properties. 

No calculus can be built without a change-of-variable formula: the Ito integrals 
provides one. Let F : R — > R be a smooth function (or at least twice continuously 
differentiable) . Also, suppose that u and v are measurable adapted processes such that 
J^u^ds < oo and \ v\ ds < oo almost surely for every t eT. For 

X{t)=Xo+ [ u{s)dB{s)+ [ v{s)ds, (3.14) 
JO Jo 

we have 

F{X{t))-F{Xo)^ I F'{X{s))u{s)dB{s)+ [ F'{X{s)) v{s) ds 
Jo Jo 

+ \jy"{X{s))u\s)ds. 



(3.15) 



This formula was obtained using the approximation by step functions we have 
previously commented. We may try guessing what happens if the point where we 
evaluate / to build the former succession, is selected in a different way. For example, 
let us take the process (3.14) and a partition 7r„ = {0 = to < < • • • < = 
interval [0,i]. The sums 



n— 1 



j=0 

converge to 



[ X{s)dB{s) + - [ u{s)ds. (3.16) 
Jo 2 Jo 

This limit integral is called Stratonovich integral. Now, comparing against the pro- 
cess X{t) itself the second term in this integral looks hke a derivative in the sense, 
^dX{s)/dB{sy. So we could write it as 



Therefore, our next question is: can such an operator be defined formally? The answer 
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is yes. It appears when one tries to define the Ito- Wiener integral for non-adapted 
processes. That is, let F{u) : 5*(M) — > R be a process such that 



^ f{{uj,(f)i),...,{uj,(f)n)), 



(3.17) 



where / e C°°(]R") has partial derivatives with polynomial growing, and the functions 
01,..., 0„ e <S(M) are fixed. Thus, we define the Prechet derivative, also known as 
Malliavin derivative, of F as 



D^F{u) = lim - + £(01, 0), . . . , w{(t)n) + £(</>„, </>)) 

£-»0 £ 

- f{w{(t)i),...,w{(t)r,))\ ; 



(3.18) 



moreover, if there exists a process DtF such that D^F = {D F, 0) — where (■, ■) is again 
the inner product in L^(]R) or L^(]R'') — we say it is differentiable. For / = id : M ^ M 
is 



0i(s) dB{s) 



(01,0), 



and thus is Dtw{4>i) = 0i(t). In general, the derivative is just the expression: 



i=l 



(3.19) 



This operator is closed and unbounded with values in L^(]R x fl) defined on the 
(dense) set D^'^ of smooth random variables with norm. 



|F||^_2 = E[|F|^] +E 



\D.F\ 



L2(R) 



contained in L'^{Q). We define the adjoint operator 6 as an unbounded operator on 
L^(M X fi) with values in L^(r2) such that: 

i) Its domain, denoted by Dom5, is the set of processes X e L^(]R x Q) with 



E 



DtFX{t)dt 



< c||i^||i,2. 



for all F e D^'^, where c is some constant depending on X. 
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II 



',) If X belongs to DomS, then 5{X) is the element of L^(n) characterized by 



E[FS{X)] = E 



DtFX{t)dt 



(3.20) 



This operator is called Skorohod stochastic integral of the process X. It transforms 
square integrable processes into random variables. It is usually written as 



d{X) := I X{t)5B{t). 



(3.21) 



This stochastic integral does not require adaptness for X; nevertheless, if it is adapted 
then it coincides with the Ito integral. Moreover, The Skorohod integral is the right 
tool to understand stochastic integrals defined by Riemann sums. 

Again, let us assume our parameter space is T = [0, 1]. It is denoted by Li^2 the 
class of processes X G L^{T x Q) such that X{t) e Di^2 for all t, and there exists a 
measurable version of the two-parameter process DgXit) satisfying 



E 



I [ {DsX{t) fdsdt 



< oo. 



This space is a Hilbert space with norm ||X||^ 2 = ll"^lli2(TxQ) 
follows that Li^2 C Dom S. 

Now, for any process X e L^(T x Q) and any partition tt = {to = < ti < • • • < 
tn-i < tn = 1} the step process 

converges to the process X in the norm of the space L^(r x fl) as |7r| = maxj \ti+i — ti\ 
tends to zero. Furthermore, it also holds in Li^2 whenever X e Li^2- This means that 
the derivatives 



n-l 



D,X^{t)^Yl 



'i+l 



U+i 



dsX{s) ds X{u,ti+i]{t) 



DsX{t). 
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On the other hand, the Riemann sum associated to the preceding approximation is: 

Thus, for any X e Li^2 we find 

(^(X'^) = - V /' /' DsX{t)dsdt] (3.22) 

moreover, it converges in L^(Q) to 5{X). Besides, this convergence does not guarantee 
the existence of the Riemann sum. Some conditions should be introduced to make the 
second term at the right-hand side converge. This summand is, in fact, an approxi- 
mation of the trace of the kernel DgXt in T^. It is undefined for an arbitrary square 
integrable kernel. The set of functions where it exists has two properties: the mappings 
s — > Dt\/sX{t A s) and s — > Dtf^sX{t V s) are uniformly continuos with respect to and 
sup5^^E[|£)sX(t)|^] < oo. Then, we have the following limits (uniformly in t): 

D+X{t) ^\imDtX{t + e) 
D-X(t) =\imDtX(t-e), 

£\0 

from it we construct the operator V = -\- D" . With all these conditions at hand 
the Riemann sum converges to the Stratonovich integral and we have 

f X{t)odB{t)= [ X{t)SB{t) + }- I {VX)tdt. (3.23) 

Henceforth, we have accomplished a definition for the Riemann 'like' approximating 
series, they are coherent with our previous views and our rough idea of derivative — see 
equation (3.16). 

3.2 Wiener-Ito Chaos Expansion and Wick product 

The chaos expansions allow us to write any given random variable as a series of 
smoothed white noise functional. There are two versions: one based on terms of Her- 
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mite polynomials, the other using multiple Ito integrals. Both version are, of course, 
related and eventually lead to the definition of a new product: the Wick product. These 
three concepts are very important, for they provide a set of analytic tools — Ito formula 
included — that will allow us to solve stochastic differential equations. 

3.2.1 Chaos expansion in terms of Hermite polynomials 

The Hermite polynomials Hn{x) are defined 

Hr^ix) = {-l)-e^y^£-{e-^y^y, n = 0, 1, • • • . (3.24) 
The first polynomials are then: 

Hq{x) = 1, Hi{x) = X, H2{x) = — 1, H^{x) = x^ — 3a;, etc.. 

Now, we define the Hermite functions — a detailed description of their properties can 
be found in Sundaram (1993): 

U^) = {2^-\n - l)\V^Y''^ e-^"/^H^_,{x)- n = 1, 2, • • • . (3.25) 

These functions belongs to <S(]R); moreover, they constitute an orthonormal basis for 
L^(]R). We will use both, the Hermite polynomials and functions, to define a basis for 

Let 5 — {5i,...,5d) e N'' denote d-dimensional multi-indices, then the family of 
tensor products 

6 — C(5i,...,5d) = 6i • • • 6d (3.26) 

is an orthonormal basis for L^(M'^). Let (^^-'^ represent a given fixed order for the set of 
multi-indices, such that, 

that is, an increasing order. Now we can define 



Vj — ^50) = Lu) Lu) ; J = 1, 2, • 

1 d 
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We will consider, in particular, the set J of all sequences a with only finitely many 
CKj 7^ 0. Therefore, for a e J" 

oo 

7^«(a;)=J]i/«,((a;,r/,)); a; G (3.27) 

These family of functions constitutes an orthogonal basis for L^{n), and 

Il^a||i2(^) = a\ := ailaa! ■ ■ ■ . 

Now, we are in conditions to formulate the Wiener-Ito chaos expansion theorem: 
every / e L^ilJ^) has a unique representation 

f{u;) ^^CaHaiuj), where c^eR. (3.28) 
Moreover, we have the isometry 

II/IIW) = E«'^- (3-29) 

Let us consider the 1-dimensional smoothed white noise as it was defined in (3.5), 
it is 

oo 

oo oo 

^J2{<P,rjj){u;,rjj) = ^(0, r?,)H,, (a;), (3.30) 

j=l 3=1 

where ej — (0, 0, . . . , 1, . . . ) with 1 on the entry number j, and otherwise. This 
convergence is in -L^(//). In this case, it is rj^j{t) — ^j{t). Also, we can calculate the ex- 
pansion for the 1-dimensional (1-parameter) Brownian motion defined in the preceding 
section. The expansion of the step function Xlo,t] , using the Hermite functions, is: 
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so the expansion for the Brownian motion is 



CO y „t s OO 

3.2.2 Chaos expansion in terms of Ito integrals 



(3.31) 



The latter expansion is equivalent to another one built using iterated Ito integrals. 
This is defined as follows: Let . . . be a symmetric function then its n-tuple 
ltd integral for n > 1 is 



/ / ■•■ / ^ti,t2,...,tn)dB(ti)dB(t2)--- dB{tn), (3.32) 

■OO J —OO J —OO J —OO 

each integrand in the iteration is adapted because of the integration limits of the 
preceding integrand. Using the Ito isometry n times whenever $ e L^(R") we find 



E 



^dB 



n\ / (^{ti,...,tnfdti---dtn = n\\\(^\ 



(3.33) 



Now, let a = (ai, . . . , q;„) be a multi-index such that n = |q;|. In 1951 Ito found a 
fundamental result: 



(3.34) 



where <S> is the symmetrized tensor product, i.e., for f,g : 



it is 



{f ®9){x,y) = f{x)g{y) 



and 



{f®9){x, y) = i^[f ®9 + 9® f]{x,y)-, {x, y) e R', 
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(the same apphes to higher dimensions). Therefore, comparing equation (3.34) with 
definition (3.27) we have 

^ e*°dS®" = 7i„(a;), (3.35) 



here we have introduced the multi-index notation = (8) • • • Cf'"*' ■ K we 

assume now that / e L'^i/j) has the chaos expansion (3.28); thus, we may rewrite / 
using the latter equation as 



oo „ 



n=0 \a\=n 

Henceforth, 

oo „ 

n=0 -^K" \a\=n 

where L^(]R"') denotes the symmetric functions in L^(]R"). Moreover, the isometry 
relation reads 

oo 
n=0 

3.2.3 The Wick product 

The representation of stochastic processes by means of the chaos expansion repre- 
sentation provides a favorable setting to study stochastic differential equations. Un- 
til now, we have characterized these processes with function and distribution spaces, 
S{W^) C L'^{W^) C S*{W^), but we will need to extend them a bit more. 

Again we impose a fixed order for the multi-index family S — {Si,. . . ,5d) e N^. 
Let us introduce the following notation: for a — (cci, . . . , a, . . . ) e and /3 = 
(/3i, . . . , . . . ) e a finite sequence it is 

= af • • • • • • where q;° = 1. 

It can be proven (Reed and Simon, 1980) that: 

i) For (f) e L^(M'^), such that = YlJLi (^jVj: where the aj = (0, 7]j) are the Fourier 
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coefficients with respect to the multi-index Hermite functions. We have e 
S{R'^) if and only if 

oo 

for all (i- dimensional multi- indices 7 = (71, ... , 7^). 

ii) Also, the space »S*(]R'^) can be identified with the space of all formal expansions 

00 

® = ^^^^ 

such that 

00 

J2b%S^'Y'' < 00 
i=i 

for some d-dimensional multi-index 7' = (7^, . . . , 7^). 

Similarly, we can define an analogue for the probability space the Kondratiev 

spaces. We will not give the more general version of these spaces, because it is not 
required in the present discussion. Therefore, let us define the quantity 

(2Nr-n(2jp, (3.37) 
j 

where 7 = (71, . . . , 7^, . . . ) e has finite non-zero numbers. The stochastic test 
function spaces Sp {0 < p < 1 fixed) are the set of all the sums 

a 

such that 

ll/IIJ := J2 4(q;!)^+^(2N)'=° < 00 for all keN. (3.39) 

a 

On the other hand, the stochastic distribution spaces S-p consist of all formal ex- 
pansions 

F = ^ baHa with bc^eR (3.40) 

a 
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such that 

\\F\\_p ■= J2 bl{a\y-''{2n)-'i°' < oo for some qeN. (3.41) 

a 

The seminorms || • ||p gives a topology for Sp, and the space S-p can be thought to be 
the dual of the stochastic test function space by means of the inner product 

{FJ) = "^baCaal. 

a 

Note that for p e [0, 1] we have 

SiGSpGSoC L^in) C «S_o C S^p C 

In particular if both F and G belong to L'^{n), then (F, G) = K[FG]. The spaces »So 
and »S_o are called Hida spaces, and respectively denoted S and S*. 
Now, we can define the Wick product: for two elements 

F — ^2 (^a'Ha, G = ^2 boHa G 
a a 

we have 

FoG^J^aabpHa+p. (3.42) 

The product is independent of the base elements of L^{ld). Moreover, the spaces Si, 5_i 
and S, 5* are closed under the Wick product. In the sense F,GeA^FoG^A with 
A anyone of the former spaces. Of course, the three laws for products — associability, 
commutativity, and distributiveness — are obeyed. 

The Wick powers F**^; k = 0, 1,2, ••• of F e 5_i are defined inductively as 
follows: 

poo ^ ^ 

pok ^ p ^ po{k-i) for A; = 1,2,--- . 

Moreover, given a polynomial p{x) — '^^=0 ^^-"^^ straightforward to define its Wick 
version, 

N 

p%F) = J2<^kF''- 

k=0 
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It can be proven that for F,G e L'^i/j) Gaussians, that is, 

F{uj) = ao + ^ Ofc^efe ('^) and G{u) = &o + ^ ^fc^e^ ('^) 

k=l k=l 

with Y.k=i E^i &fe < oo, it is 

oo 

(G'oF)(u;) = {GF){u)-Y,(^kbk. (3.43) 

fc=i 

Where it had been used the property 



Applying this formula to the smooth white noise expansion (3.30) we find 

w{(p) ow{ilj) = w{(p)w{i') - {(f),tlj); (3.44) 
moreover, ii t/j — (f) and \\(j)\\ — 1, then we have w{(f)y^ — H2{w{(j))), and in general: 

w{<f>y^^Hr,{w{<f>)). (3.45) 

3.2.4 Skorohod integration and Wick product 

The Skorohod integral can be written in terms of the chaos expansion. Let Y{t) = 
Y{t, uj) be a stochastic process such that E[y (t)^] < oo for all t. We already know that 
this process can be expanded as 

oo „ 
n=0 -^^^ 
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where /„(•, t) e L^(M") for n = 0, 1, 2, . . . and for each t. We denote by . . . , 

the symmetrization with respect to the n + 1 variables. Thus, assume that 



L2(R"+i) < OO. 

n=l 

We can define the Skorohod integral of Y{t) as 



(3.46) 



It has the norm 



Y{t) 5B{t) 



= J^(n+ l)!||/„||L2(Kn+l). 

L2(K) „=o 



On the other hand, we say Z{t) — Ca{t)Ha £ S* is <S*-integrable if from its 
chaos expansion the expression 

/ Z{t) dt^Y^i I Ca{t) dt] Ha{0j) 

belongs to S* . Now, the process 

OO 

fe=0 

because the Hermite functions are bounded: ^n(^) < rT^^^"^ . From equation (3.31) we 
have ^ 



Therefore, we have proven that dB(t)/dt = W{t) is well defined in S*. Afterwards, 
a last fundamental theorem remains to be addressed (Holden et al., 1996, Theorem 
2.5.9): assume that Y{t) = ^^^cJi-ci is a Skorohod integrable process, and let a < 6 
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real numbers. Then Y{t) oW{t) is <S*-integrable and 

b pb 

Y{t) 5B{t) = / Y{t) o W{t) dt. (3.48) 

J a 

3.3 Stochastic Calculus for fractional Brownian mo- 
tions 

In the past years different approaches have been given to produce a Stochastic Calculus 
for fractional Brownian motions: Zahle (1993,2001), Decreusefond and Ustiinel (1998), 
and FoUmer et al. (1995). Basically, these approaches tackle the problem of construct- 
ing a calculus, but from two different starting points: one uses a pathwise definition 
of the integral while the other rests on the Malliavin Calculus as we sketched earlier 
in this chapter. In all these circumstances the processes have are persistent. We will 
follow Hu and 0ksendal (1999) and Duncan et al. (2000) into the second approach. 
We will construct a Stochastic Analysis from a Chaos expansion. 
Let : R+ X R+ ^ R be defined as follows 

0(s, z) = H{2H - 1) |s - zf^~^ , 

for a fixed H E {1/2,1). Then we say / e -f'^(R) if it is measurable and 

\f\l ■=11 ns)f{z)4>{s,z) dsdz < oo. (3.49) 

Afterwards, the inner product can be defined in L^(M), 

if,9)<t>-= [ [ f{s)g{z)cf>{s,z)dsdz, forall/,^ eL2(M); 

therefore, I/^(M) becomes a separable Hilbert space. 

Again, we take 5(M) C L'^i^) to be the Schwartz space of rapidly decreasing 
smooth functions on M. Its dual fl = S'{M.) is the probability space with the associated 



J n. 
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probability measure, /z^, found applying the Bochner-Minlos theorem, 

Jn 

where {cu, f) is the usual pairing between elements in the dual and functions on R. 
Because of the latter construction this probability measure can be shown to induce 
properties like those in (3.7), i.e., 

E[(-,/)] = 0, and E[(-,/)2] =1/1^. (3.50) 

Once more, the triplet (f2, B{fl) , /x^) becomes a probability space — B{fl) is the Borel 
algebra on f2. It is usually called fractional white noise probability space. 

Now, let L'^{iJi(j)) — L'^{Vt,B{Vt) ,//^) be the space of all the random variables X : 
n ^ K such that 

ll^lli^W):=E|X|2<oo. (3.51) 

Furthermore, the functions in I/^(M) define a set of random variables of the form 
/(a;) = (cu,/). It is included in L^(/x^); that is, the condition (3.49) induces square 
measurable random variables because of equations (3.50). 

With the same arguments as before we have that: 5(R) is dense in -£/^(IR); for any 
/ e (M) the series /„ G M are such that /„ ^ / in -^^(M); and so, the following limit 

lim(o;,/„) (3.52) 

n— >oo 

exists in L^(/i<^). 

We define now the fractional Brownian motion process as follows: 

B"{z) :- B^'iz^u) = (a;,X[o,.)) e L^fi^). (3.53) 
For simplicity we will thought designates the 2;-continous version of the rightmost 
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hand side term. As for the step function Xlo,z) : M — >^ [— 1, 1] again: 



XlO,z)( 



1 if < s < z 
-1 if ^ < s < . 

otherwise 



Because of property (3.52) we have again that for any / G i>|(IR) definition (3.53) is 
equivalent to 

{cuj)= [ f{z)dB\z,uj). (3.54) 



Under the same procedure we can verify for f,g& that 

n{u,f){u,g)]^{f,g)^. (3.55) 
For / as above, we define the exponential function £ : I/^(K) L'^{ii(j)) as 

£(/) = exp^/ciS^-i|/|^). (3.56) 

Thus, the Hilbert space i^|(M) is naturally associated with the fBm process from the 
formulation as an abstract Wiener space. Let 8 be the linear span of the exponentials, 
i.e., 

£ = \^au£{fk)]n e N, afe e R, h e Lj(R) for A; e {1, . . . , n}| , (3.57) 
is dense in L'^{ii^). 

Nevertheless, some tools we are going to introduce here require a more familiar 
functional expansion, and the Hermite functions (3.25) will help us again. First, we 
note that we can map the orthonormal basis they form in L^(R) to an orthonormal 
one in -L|(R) through the isometry map (see Lemma 2.1 in Hu and 0ksendal, 1999) 
L = r^^^„ defined 

{z-s f-'/^fiz) dz, 
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where 

'H(2H - i)r(f - H) 



ch 



r{H - i) r(2 - 2H) ■ 

Prom the identity (Gripenberg and Norros, 1996, p. 404) 



/zAs 
-oo 



we see that 



z 



C„(S) Z) ds = CH I {Z- Sf~^'^in{s) ds 



— because the ^„'s are an orthonormal basis these integrals are also smooth. 

Let X be the set of all finite multi-indices a — (ai, • • • , a^) of nonnegative integers, 
we define 

In particular, if we put a = e* then, in the very same way as in Section 3.2, we get 
from (3.54) and the definition of Hermite polynomials 

HA^) = Hi{'^,ii)) = = I Us) dB^'is). 

Jr 

These functionals are elements of LF'{^^)^ and they form its basis (Duncan et al., 2000, 
Theorem 6.9). That is, for X G L'^{lJ'^) there are G M and a G X, such that 

X(u;) = J]c„7i„(a;), (3.58) 

and also 

m\w,) = Y.^'-^l- (3-59) 

These coefficients are given by Cq, = E[X?ia] 

The existence of this property let us define a fiavor of (fractional) Hida spaces: the 



3.3 Stochastic Calculus for fractional Brownian motions 



102 



fractional Hida test function space Sh which is the set of all 

ip{u!) — ttaTi-ai^^) £ L'^{l^(j>): such that 
1^^^ = ^Q;!a^(2N)'=" < oo, for all A; eN, 



where 

(2 = \{{23y^ for any element 7 = (71, • • • , 7^) e X; 
3 

and the fractional Hida distribution space S^, the set of all formal expansions 
Y{(jj) = ^^6^?i^(a;), such that 



(3.60) 



Using these definitions is not hard to see that Sh C L'^^jJ^d}) C S*^. 

It is now time to show how the fractional white noise and integration with respect to 
is defined. Let us first calculate the expansion for the stochastic integral in (3.54). 
For any / e -^^(K) — any given deterministic function — we have from equations (3.58) 
and (3.55): 

« 00 

/ f{s) dB^'is) = 5](/,|,)0K^(^)- (3.61) 
When / = X[o,z) in the left hand side we recover (3.53) and the following relation holds 



B'^iz) -Y.\f [I Us)(t>{s,u) ds\ du 



(3.62) 



if we check its norm 



\B"{z)\\l-, = E / ( / ds) du 



(2A;)-« < 



k=l ^ ^ 



(3.63) 



3.3 Stochastic Calculus for fractional Brownian motions 



103 



(C is the Riemann's zeta function) because 



s)4>{s, u) ds 



ch 



J (— oo,u] 



<Mk^/^, (3.64) 



here we have used the bound for the Hermite functions given by Szego (1967, pp. 
198-201). Furthermore, when q > 4/3 the former inequahty also shows that is 



continuos and differentiable in S^. Its derivative 



d 
dz 



B'^iz) = y^{l ik{s)4>{s,z)ds\ n,.{u) := W'^iz) e S*h (3.65) 
k=i ^-^^ / 



is the formal definition of fractional white noise. This noise is also continous in S}j, 
when z > s 



W'iz) - W''{s)\\%_^ = Y,e'\ / Uu)(t>{z,u)du- / Uu)(t>{s,u)du 

< E i [(^ Ms + u) I du 

k=i L-^o J 



{2k)-'' < 



(2A;)-« < 



22-«4m2 / 1 , , 



(3.66) 



the same holds when z < s. 

Of course, this chaos expansion has its own Wick product. Let X{u!) — X^^gx ^a^ai'^) 
and Y{uj) — Ylip^^fi^pi'^) be in S*jj^ then 



{XoY){uj) = aahpHa+pii^) = E I E ""^^ ) '^li^)- 

a,i3ei -yei \a+/3=7 / 



(3.67) 
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For example, let f,g & L'^i^) then using equation (3.61) we find 



J=l / \j=l 

oo 



i,i=l i=l 

' oo \ / oo \ oo 

I f{s)dB"{s)\ ■ ( I 9{s)dB"{s)\ - {f,g)^. (3.68) 

This property is a special case of a more general one for Gaussian variables, that is, 
for X ^ ao + J2'i^o (^i'Hei and Y ^bo + J2^o hjHei we have X oY ^ X - Y - Y^^^i (^tk 
as was proved in (3.43) for the Brownian case. Afterwards, ior f — g — proceeding 
inductively with the latter equation yields 

Now, as we extended polynomials into the Hida space for Brownian motions, we do 
the same here but with the power series. The Wick exponential defined by the power 
series 

00 

expO(X) = -r^'", 

n=0 

provided it converges in S^. It has the same algebraic properties as the usual expo- 
nential, e.g.: 

exp*(X) o exp*(F) = exp*(X + Y). 
This Wick exponential is the keystone of this section, for it provides a hnk between 
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the two expansions given here. If we set X — a{u!, ^j), it is 

^-^ nl 

1=1 

°° n 
i=l 

= exp (^a{uj, ii) - ^a^^ , (3.69) 

because of 

/ 1 \ ~ 
exp \tx--t 1 =2^^Hn{x). 

^ ' i=l 

Therefore, when X = {u!,f) = J^fdB^ 

/ oo 

exp*((u;,/)) = exp* I ^{f,ii)4,{uj,l 
\i=i 

o 
i=l 

°° / ~ ~ 1 ~ 



i=l 



(oo ^ oo \ 

i=l i=l / 



and thus 

exp*((..,/))=^(/), (3.70) 

as the right-hand side was defined in (3.56) the relation between the expansions is 
settled. 

It is appropriate to show now the behavior of the Wick product within an average. 
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Let X = ^j^gj aaTi-a and Y = JZp&i bfiUp have the usual chaos expansion thus 

7eX \a+/3=7 / 
7eX \a+/3=7 / 

= aa6^0! = ao6o = E[X]E[y], (3.71) 

a+yS=7 

here we used the fact that the TY's are an orthonormal basis. 

Obviously the next step is to introduce the (fractional) MaUiavin derivative for 
these processes or 0-derivative, for X e L'^ifJ-^) and g G -t'^(IR), the alternative version 
to (3.18) reads 

D^gX{u;) = liml\x{u + 6 [ {<^g){u) du)) - X{u)\ , 

where {^g){z) = J^(j){z,u)g{u) du. Afterwards, if there exists a function D'^{s)X such 
that 

D^gX = / {DfX) g{s) ds, Wg e (M), (3.72) 

we say that X is 0-differentiable, and DfX is the 0-differential. Let us point out some 
properties for the fractional MaUiavin derivative, with X defined be as always and 
/, : R — > R these are: 

D^JiX) = f{X)D^gX, (3.73) 
D^g{u;,f)^{f,g)^, (3.74) 

Dt{uj,f)= [ 4>{u,s)f{u)du. (3.75) 
Let us inspect another property for this operator. We can compute the second 
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moment of S{f)S{g). Because of IE[^(/) <>^ig)] — 1, 

1 = E[s{mg)] - J2 ^,m^,f){oo,9)r] 

n. 

n=l 

n. 

n=l 
n=l 

^E[S{f)S{g)]-eMf,9)<t>+h 



we used property (3.55) in the last steps; so, E[S{f)S{g)] = exp(/, g()<^. We construct 
the following, 

E[{£{h) o 8{5f)){£{h') o £{6g))] = E[£ih + 5f)£ih' + eg)] 

= exp(/i + 5f, h' + eg)^. 

Taking partial derivatives in 5 and e afterwards yields 



E 



8{h)o f fdBA(s{h')o f gdB" 

= exp(/i, /i% [(h, f)4,(h', g)^ + (/, g)^\ 

= E[D^^S{h)D^^S{h')+S{h)S{h'){f,g)^] 

Henceforth, because any two X,Y & L'^il^cp) can be decomposed by the span £ we 
finally find. 



EU Xo 



I mdB^is)) (yoI g{s)dB^{s) 

= E[{D^fX){D^,Y)+XY{f,gU . (3.76) 

This equality will allow to change the integrator inside (3.54) by a stochastic function 
X : R X Q — > R such that < oo. That is, define the stochastic integral for 

fractional Brownian motion. 

The basic procedure consists of building a Riemann sum, replacing the standard 
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n-l 



product by the Wick one, 

1=0 

Observe that for any partition tt = (2:0 < zi < ■ ■ ■ < Zn-i}, 



(3.77) 



E 



n-l 



J2x{z,)o{B^{z,+^)-B^{z,)) 



i=0 



n-l 



^J2^[X{z,)o{B^{z,+,)-B"{z,))] 

i=0 

= j2m{zi)]mB"{z,^i) - B"{z,))] = 0. 



i=0 

Next, we compute the L?'{^^) norm of the former sum. Note that, 

E{ [X{z.) o {B^{z,^,) - B"{z,))] [X{z,) o {B^{z,^i) - B^{z,))] } 

/ DfF{zi)ds Df;F{zj) du + X{zi)X{zj) / (f){s,u) duds 

J Zi J Zj J Zi J Zj 

is obtained from (3.76); afterwards, 

^[{Sn{X)f] = Y.A ^tF{zr) ds / DtF{z,) du + 
i,j=o Zi J Zj 

+ X(zi)X(zj) I I (f){s,u) duds 

The continuity of X and the existence of the trace of DfX (Duncan et al., 2000, 
Theorem 3.9) makes this sequence converge in -^^(/x^) as |7r| — > 0, and it converges 
to 

|2 



E 



^ DfX{s) ds^ +\X\ 



In these conditions we say it is the fractional Brownian Stochastic Integral: 

-L 



hm SniX) ■= [ X{s)dB"{s)] 

n-»oo Jo 



(3.78) 



3.3 Stochastic Calculus for fractional Brownian motions 



109 



moreover, the following equality holds 

[ X{s)dB"{s)= I X{s)oW^{s)ds, (3.79) 
JO Jo 

while the integral on the left-hand side represents the hmit (3.78), the right-hand side 
is just the integral evaluated under the Hida expansion of the Wick product defined in 
(3.60)-(3.67). 

Dropping the Wick product in definition (3.77) still produces a limit if the conditions 
given above are satisfied. This integral is the Stratonovich integral X[s) o dB^[s) 
, because 

n-l 
1=0 

n— 1 n—1 

= J2 X{z,) o (i^^(^.+i) - B'^iz,)) + ^ D^^^^^^^^^^^X{z,) 

i=0 i=0 

n-l n-l ^^^^ 

= J2xizi)o{B^{z,+i)-B^{z,)) + J2 / dsDfXiz,), 

we have 

pL pL pL 

/ X{s) o dB"{s) = / X{s) dB"{s) + / D'^^X{s) ds. (3.80) 
./o Jo Jo 

This property is the counterpart from (3.23) in the Brownian motion it should 

be if the analogy follows from (3.48) into (3.79). But here both operators on the right- 
hand side can be evaluated without difficulty. We will finish this chapter with three 
theorems from Duncan et al. (2000) we will employ soon: 

Theorem 4.2 LetX{z) a stochastic process defined as above, anc? sup^g [q E |Z}|'X|^ < 
oo. Also, letr]{z) = X{s) dB^{s). Then for s,z e [0,L) 

Dfri{z)= / DfX{u)dB"{u)+ / X{u)(j){s,u) du. (3.81) 
Jo Jo 

Corollary 4.4 Let r], = f{s) dB^{s) and F{z, x) -.R+xR^R, where / e L; 
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is continous and F has second continous derivatives. Then 

^{s,n{s))ds + ^—{s,rj{z))f{s)dB"{s) 
+ ^{s, vis)) Hs, s')f{s') ds'ds. 



(3.82) 



Theorem 6.11 If X & L'^{fJ'(j>) then there exists a sequence {/„ e L^(]R")}„gN such 
^hat Yln=i \fn\l<oo and 



oo „ 

X = E[X] + V / Us^,--- ,Sn)dB^^---dB^^ 

n=l M 



(3.83) 



where 



/ fn{si, ■■■ , Sn)fn{s'i, ■■■ , s'J 0(si, s[) ■ ■ ■ 4){Sn, 4) dsi'-- dSn ds[ ■ ■ ■ ds'^ 



L^(R") is the n-dimensional space of symmetric functions. Given the base complete 
orthonormal base {^n}neN C L^(]R+) then L|(R") is the completion of all function of 
the following form: 

/(si,...,s„) = ^ afei,..,it„4i(si) 42(^2) •••|fc„(sn) ■ 

l<fcl,...,fen<fc 

Associated to this functions we define the multiple integral as 

l<ki,...,kn<k "^^ 

It is not difficult to prove (Duncan et al., 2000, Lemma 6.6) that given / e L^(M") 
and / e L2(M™) it is 

. if, Q)d> it n — m 
nW)U9)] = { ' ' , , ■ (3.84) 

if n 7^ 777, 
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Moreover, for the iterated integral 



J 0<Sl<S2<---<S„<t 



f (f Us,,--- ,s^_,,s^)dB^^...dBl_\dBl 

Jo \J0<Sl<S2<-<S„-l<Sn J 



is n\ times In{f)- 



Chapter 4 

Stochastic Geometric Optics 



Diverse experimental techniques have been devoted to the study of the optical prop- 
erties of the turbulent atmosphere. Plenty of them are based on the analysis of the 
output of laser beams making their way through it. But also, controlled experiences 
had been developed for the laboratory, such as the experiments performed by Con- 
sortini and O'Donnell (1991, 1993); Consortini et al. (1997). These experiences apply 
Geometric Optics to interpret the data acquired. 

All these studies have their theoretical grounds on the precursor paper by Beckman 
(1965), who was able to find a nice relationship between the variance of the turbulent 
refractive index fi{r) — being homogeneous and isotropic — and the variance of the laser 
beam wandering over a screen. As it was pointed out in Chapter 1, he proposes (1.100) 
as covariance function because it gives meaning to the derivatives of the refractive 
index. Moreover, he pointed out that the Kolmogorov-like structure functions ". . . are 
mathematically fairly unmanageable". The literature after him forgot this warning: 
modifications to his solution were given (e.g., Consortini and O'Donnell, 1991) but for 
the wrong covariance, the Kolmogorov structure function. 

We do intent to show here, armed with our refractive index's model, that the ray- 
path equations are manageable. But this requires the Stochastic Calculus we have 
introduced in the last chapter. 
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4.1 Introduction 

Before start working in our approach we will briefly describe the differences between 
it and other works. In most of them the markovian model plays a central role. Le- 
land (1989) exaustively depicted by it. The markovian model provides the following 
covariance: 

E[ 6(p; ^)e(p'; z')] = 5{z - z')A{p - p'), (4.1) 

where A is a differentiable function as defined in Appendix B. This covariance is asso- 
ciated to a process build from the Brownian motions' distribution space to a bounded 
linear operator L on some Hilbert space H] that is, e = L{B^/^). Since this opera- 
tor can be described by using some kernel function whose coeficients are differentiable 
functions in p. Obviously, this model transfers all the discontinuities to the z-axis. 

For instance, let us illustrate the problem with the simple example: choose L{B^/'^) = 
Jq F{p; s)5^/^(s) ds. Assuming F is continuously differentiable in p, the following 



F,{p-s)dB'l\s) 



(4.2) 



is well-defined. On the other hand, the covariance of the original process is 



E 



'l{B^''){P-z)L{B^/'){p'-z') 



zl\z' 



F{p,s)F*{p',s)ds; 



therefore, differentiating the above by , we find 



dxdx' 



dxdx' 



■E 



L{B'/')ip;z)LiB'/'){p';z'] 



zAz' 



F,{P,s)F;{p',s) ds 



Henceforth, from equation (4.2) we observe that 



a^E 



L{B'l^){p-z)L{By''){p'-.z') 



dxdx' 



E 



^L(5V^)(p;.)Al(sV2)(p'.,') 



dx 



(4.3) 



It is this property the commonest property used in turbulent optics not regarding its 
original nature; that is, equation (4.1) or the like. Moreover, we can also evaluate 
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the fractal dimension of this type of processes. Let us use the Kolmogorv's criterion^ 
(Kunita, 1997, Theorem 1.4.1, pg 31.) for that. Let n be an even integer, it is 



E< 



Jo Jo 



E 



Jo Jz' 



Now, if wc name Gs = G{p,p',s) = F{p,s) - F{p',s) and Hg = F{p' , s)x[z',oo){s) 
after applying the Newton's binomial theorem, then we will have a summatory with 
the following terms 



^^E[/,(G)4_,(F)], with Ir,U)=(^f^fsdBll'^ , 



these integrals can be turned into symmetric integrals as the ones shown in the latter 
chapter. We note from the orthonormal property of stochastic symmetric integrals — an 
equivalent to (3.84) for the Browinian case — that the only remaining are three: 



E[4(G)] =E[/^/,(G)] = {n/2)\\r {F{p,s)-F{p\s)f ds 

Uo 

Wn{F)] = E[/^/,(F)] = (n/2)! [ f F\p', s) ds 

-J z' 



E[4/2(G')4/2(F)] ={n/2)\ 



r{F{p,s)-F{p\s))F{p',s)ds 

.'J -2' 



n/2 



-'^ Given a process ^(r), with r in a closed domain D in R''. Assume that there exist positive 
constants s, M and a^, i = 1, . . . , with o.Q^d = J2i=i < ^ satisfying 

d 

E[|X(r) - X(r)n <M^\xi- x'p , for every r,r' e D. 

i=l 

Then it has a continuous modification X such that 



X(r) - X{r) < K{lo) ^ \x^ - x'f' , for every r, r' e [0, l]'', 



holds for almost all w. The coefficients /3i are arbitrary positive numbers less than ai{ao — d)/aos. 
We call it a . . . ,/3d)-H61der continuous process. 
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It is not hard to find bounds to these, 



(n/2)! 



{Fip,s)-F{p',s)f ds 



n/2 



< (n/2)!Mi||p-p'r, 



(n/2)! 



/ F\p\ s) ds < (n/2)!M2 \z - z'f , 

Jz' 



(n/2)! 



r{F{p,s)-F{p',s))F{p',s) ds 

Jz' 



n/2 



< {n/2)\M4p-p'f/^\z-z' 



I n/2 



Finally, using the property \\p — < 2"/^(|x — x'l" + |y — y'|"') we have 



{ [l(sV2)(p; z) - L(sV2)(p'; ^')] "} < C (|x - + b - + N - ^T^') ■ 



(4.4) 

Therefore, we observe that /3i^2 < (n — 4)/n and Ps < |(n — 4)/n. In particular, 
min{/3i, P2, P3} — P3 < i(^ ~ 4)/n < 1/2. Using the Holder continuity we observe this 
process gives a isoscalar fractal dimension less than inf{3 — Ps} < 2|. Moreover, the 
fact that 



m^llr - r'lP < E 



r F{p,s)dBl/'- r F{p',s)dBl/' 
Jo Jo 



a 



provides us a bound for the potential theory, and thus we will obtain — as we did in our 
first chapter — a isoscalar fractal dimension equal to 2i. 

Therefore, not only this model does not match the covariancc function but also 
does not provide the right dimension for the refractive index. It effectively allows some 
degree of differentiability but at the cost of eliminating some physical informatin from 
the refractive index covariance. Moreover, this markovian approach is not isotropic, 
and an isotropic version will inexorably lead to a non-differcntiable process. 

In particular, wc may cite the work of Consortini and O'Donnell (1991). They 
follow Beckman's steps to evaluate the covariance of the displacements of a ray over a 
screen. Ending up with an equation of the form 



Ax^ A f j ^dzdz', 

Jo Jo 
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Figure 4.1: The graphic shows the behavior of the log of the variance against the 
distance L. Interpolating lines can be calculated and the values of their tangents are 
shown. 



where A is some constant. Afterwards, the authors commutate the derivatives with the 
average. But they do not mention the markovian approximation as the cause of this, 
and soon after they replace the covariance function by the isotropic one. This violates 
the valid use of the commutation property (4.3); since an isotropic process does not 
provide derivatives, the above equation has a priori no meaning thus it is not true we 
can commute operators. 

Moreover, we observe the markovian model is dependence on the characteristic 
lenght L as L^/^, then the former integral behaves as L^2. The covariance of the 
displacements will grow proportional to L^. This is a quality of the Brownian or 
markovian processes. 

Finally, we observe in Figure 4.1 several plots of the logarithm of displacement 
covariance against the distance, based in the experimental data found in the work of 
Consortini et al. (1997). In all cases the estimated power is below the theoretically 
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estimated. Just in the higher cases the error is wide enough to cover the calculated 
value a — 3 and its value near it. 

Next, we will use the isotropic fractional Brownian model within the Geometric 
Optics to obtain an equation for the rays. We will show that under the correct frame- 
work a solvable stochastic equation exists and its result can be directly applied to the 
problem of a ray wandering over a screen. 



4.2 Stochastic Differential Equations in Geometric 
Optics 

4.2.1 The ray-path equations 

As it is well-known, the Fermat's Extremal Principle is in the foundations of the Geo- 
metric Optics, that is, to find the ray trajectories we must find the variational solution 
to 

d( [ nds] =0. (4.5) 



We shall denote this solution by q{T), and r is a parameter with, in principle, no 
physical meaning. In Optics Treatises this parameter is usually replaced with one of 
the trajectory coordinates, which fulfills dqi/dr > 0, and is thus called the propagation 
direction. But the election of this parameter can not be done at will (Synge, 1937), 
since, for any parameterization chosen, the Optical Lagrangian 

L{q,q) =n(g)||g||, 

(g, g e are the position and velocity respectively^ ) is degenerated. It is easy to show 
this property. Calculating the momentum. 



is the configuration space. Usually is denoted by Q. 
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we see that the Lagrangian is rewritten as, 



(4.7) 



Since it is homogeneous in the velocities we can recalculate the momentum and find, 



for any pair (g, q): this matrix is singular. As it is proved by Marsden and Ratiu (1999, 
Theorem 7.3.3), the solution is not univocally determined because the second order 
dynamics equation 



obviously, can not be built. Nevertheless, equation (4.6) provides us more information, 
for it induces the following relation 



which indicates that the choice of coordinates and momenta is not free. 

The degeneracy of the Lagrangian should be worked out in the Hamiltonian frame- 
work because of the constraint we have just found. This problem of constrained Hamil- 
tonians is known as Dirac's problem in the literature. The procedure is to reduce it to 
a Lagrangian problem: because given a set of constraint functions 




Therefore, 






(4.8) 



^^1(5,?) = 0, . . . ,V'fe(g,g) = 0, for some {q,q)eTQ, 



associated to a Lagrangian L, there is a solution q : [a,b] —>■ Q (critical point) if and 
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only if 3Afc : [a, 6] — > R such that the following equation holds 



d /9L\ dL 
dr \dq^ J dq^ 



d / dijjj \ dijjj 
dr \ dq^ J dq^ 



^ dq' 



(4.9) 



(see Arnold ct al., 1993, for a proof). Afterwards, we apply this theorem to L{p, q, p, q) = 
^(p? q) ~ H{p, q) — where we have chosen the configuration space P = T*Q, and 9 is the 
canonical 1-form on T*Q — with 



^^{p,q) = 0,...,^k{p,q) = 0, 



(4.10) 



for (p, q) e P. We thus find, from (4.9), 

. d 



dp 

d 



(4.11) 



Also, we can calculate the dynamics equations for the constraints (4.10), that is. 



^, = {vl/„i/} + ^A,{vl/„vi>^.}, 



(4.12) 



where {•, •} is the Poisson bracket. The set of these equations is called compatibility 
condition set: if {^'j, ^j} 7^ then the multipliers Aj are uniquely defined. Otherwise, if 
some {^i, "ifj} are zero we have a new set of constraints, called secondary constraints, 
that should be added to the original constraints. But, when we have k — 1 and 
{^i, H} — then Ai is arbitrary. 

Now going back to our problem, equation (4.8) provides us with the constraint 

nP,'l)-l[\\pr-n'{q)], 

and the Hamiltonian H obtained from the original Lagrangian is, combining equations 
(4.6) and (4.7), 

H = ^Piq' -L^ ^Pig^ - ^Piq' = 0. 
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We just need to build the new Hamiltonian, as (4.11) suggests, 

H{p, q) := H{p, q) + A*(p, q) = A*(p, q). 



By doing so, we obtain the following dynamic equations 



dH 
dq 



dq 



A^ 



. dH , 
op op 



(4.13) 



and the constraint. 



= *(p,g) = - [IIpI 



|2 2 

— n 



(4.14) 



Finally, to ensure A is well defined we have to check the compatibility conditions. 
Because our original Hamiltonian is zero, {H, \1'} = 0. The constraint is arbitrary; 
moreover, it is actually a smooth function on the constrained space that can be freely 
chosen. There are no secondary constraints derived from the compatibility conditions 
so (4.13) and (4.14) completely define our problem (Blagojevic, 2001). 

Combining the pair (4.13) of Hamiltonian equations yields to the following second 
order equation: 

d f 1 dq\ A^ 2/ 



dr \XdT 



W,n\q{T)) 



with 



AV(g). 



(4.15) 



(4.16) 



We observe that with each selection we make for A the parameter r is also set, i.e. 
if we choose A = n^^ then 



and 



ds = \\q\\dT — dr 



(4.17) 



r is then the arc-length. But selecting A = 1 gives us ds — ndr and now the parameter 
is T = J ds/n. 
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4.2.2 Linearizing the trajectory equations 

The ray equations we have just found are evidently nonhnear, so in this section we are 
going to Unearize them. But first, we must define the parameter r and the refractive 
index. Let n be the refractive index of the medium and no its average, as it was defined 
in Section 1.3.1, we write 

n^{q) = nl + a€*{q), (4.18) 

we changed the stochastic permitivity €{q) by ae*, where 0{e*) ~ 1, so the strength 
of perturbation is due to a. This term contains all the inhomogeneities of the media, 
thus when a — the index is constant. Now, we suppose the solution to (4.15) can be 
expressed as power series on a, i.e.. 



oo 



?(T)=?0 + $^«"?n(T). (4.19) 

n=l 

Also, we should develop a series for the constraint function A. Instead of using an 
undetermined constraint we will set its value beforehand: from all the possible param- 
eterizations we choose the arc-length (4.17). Now, we can rewrite equation (4.15) as 
follows 



dr^ 2 

therefore, it is better to expand 



aX'V,e* + 1 (V,A2 • q) 



in short we will write := (-l)"e*"(g)/no""^^ — note that 1/A^ is exact. 

Now, we must insert both power series in a, the expansion (4.19) and the latter 
for A^, into (4.20). We will obtain afterwards a family of differential equations from 
claiming the equality between the coefficients on the right and left for the same power. 
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The second term on the right-hand side is tricky, 



A2 



I v-^ kVe* ^ 2 k 



00 n 



Tin 



yn=0 



a' 



n=l k=l 







n=l k=l 

- (Ve* ■ %) a+ 



a 



n+l 



n=l fc=l 



n=2 



n-1 



k=l 



k=l 



= - (Ve* ■ qo) a+ 

n 

J2 kXUiVe* ■ qn-k) + E(^ - • 4 



n-l 



k=l 



n=2 
n=l 

Thus, we finally have: 



• qn-k 



k=2 



a . 



and when n > 2 



d?qn 1 



<fqi_J_ 
'd^ ~2^o 



m=l 



-A 2 * dqra-k 

.fc=i ^ 



dq„ 



dr 



a 



(4.22) 



(4.23) 



With the same criteria we obtain a constraint condition from (4.16) for each differential 
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1, (4.24) 
0, (4.25) 



equation above; 

dr 
dqi dqo 
dr dr 

y ^ • %^ = 0, for all n > 2, (4.26) 
^ dr dr ~ ^ ^ 

while the first constraint normalizes the zero-order solution, the second establishes it 
is orthogonal to the first-order solution. 

We can readily find the solution for the zero-order equation in (4.22). The result is 
the linear relationship: Q'o(t) = ar-l-b. Given that the initial condition to the problem 
is 

q{0) = 0, (4.27) 

it implies that b = 0; also, using the constraint condition (4.24) we obtain ||a|p = 1, so 
we are free to choose the coordinate frame best suited to our purposes. Let us choose: 

ze^ := qo = re^, (4.28) 

this will be our forward direction of propagation. Now, we proceed to calculate the 
next differential equation: the first-order constraint condition (4.25) reads then 

^ = 0. (4.29) 

dz 

This and the initial condition (4.27) make the component along the 2;-axis null all 
over the ray trajectory. Of course, this constraint is compatible with its corresponding 
dynamical equation (4.22). Therefore, at first-order in a we just have a differential 
equation for the perpendicular (to the direction of propagation) displacements: 

d? a 

j^Q- ^ V,e* {ze, + Q + . . . ) , (4.30) 

where Q = aqi. If we want to introduce the model we have previously introduced 
we just make the change e* , and so the parameter a — 2 l^^fA^. From the 
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values given to structure constant and the inner length, in the ideal case, we estimate 
a ~ 10~^. Afterwards, in order to examine the stochastic behavior of a wandering 
beam it will be enough to consider this first-order equation. 

With the tools we have used until now further analysis can not be done: the proper- 
ties of the turbulent refractive index must be introduced in order to completely linearize 
the former equation. 

4.3 The Stochastic Volterra Equation 

As we already know, the gradient in equation (4.30) should be given when looking for 
a solution; thus, we must provide a context to understand the previous equation. That 
is, a stochastic equation is not only determined by the type of process (the fractional 
Brownian motion in our case) attached to it, but also by the integro-dijferential theory 
employed in defining its derivatives. Moreover, there are distinctive stochastic integra- 
tion methods whether H > 1/2 or H < 1/2 (Decreusefond, 2000). Here we are going 
to make use of the stochastic calculus exposed in the last chapter, so only the if > 1/2 
case will be considered. By doing so, either we are considering the inertial-diffusive 
range, in the following sense 

^ 1 
U> 3, 

or the anisotropic scalar situation Cn — 1 (Elperin et al., 1996). The physical interest 
about this particular situation comes from the many optical experiments where as- 
pects regarding the creation of turbulence are neglected. Usually, heaters are used to 
create a turbulent medium but neither buoyancy or the temperature distribution are 
measured nor controlled, opposed to the conditions we have given through this work. 
Furthermore, the isotropic state of the index can be questioned. 

Afterwards, because the turbulent refractive index oscillates around its mean value, 
it is expected that the light wanders around the z-axis over the screen (which corre- 
sponds to the case a = 0). So, the solution we are looking for should also have 
expectation zero. This is easily achieved by the formalism we introduced: the frac- 
tional Ito integrals have expectation zero as it is seen from properties (3.53), (3.71) 
and their definition (3.77). Henceforth, using the model's definition (1.103) we can 
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calculate the gradient of the refractive index: 



[B»(l, 



E 

oo 



JL 



^k{s)(p{s,u) ds 



ik{s)<p{s, u) ds J du 
d 



dxi 



(/oir||)7i,.H 



x\ 



for i = 1,2. This equation should be understood within S"^ — it has nothing to do with 
the usual concept of derivative: we have used the chain rule and the fractional white 
noise definition (3.65). 

The procedure to interpret equation (4.30) requires to replace all the ordinary 
products containing stochastic variables by Wick products. If we do not follow this 
rule, the integrals should be interpreted as Stratonovich integrals. Thus, we observe 
from (3.80) that the mean value of the solution is non-zero, and we do not want that. 
Henceforth, 



dz^ 2 IqUq 



W^[l^^\\ze, + Q\\) 



Q\\ 



oQ. 



(4.31) 



Still, besides the changes, we have a non-linear stochastic differential equation. Worse 
than that, we have a composition of two stochastic processes. We have to find a 
reasonable way to define it. In the last chapter we explained that because any analytic 
function is expressed by a power series, it can be extended into the Hida space — 
whenever a stochastic process is an argument for it — by replacing the powers by Wick 
powers. We are going to extend this substitution rule. The representation for the 
noise in is a series with analytic functions as components (3.65); thus, it is valid 
(0ksendal, 2002, private communication) 



/ (l){s,z)ik{s)ds ^ / (f)''{s,Z)ik{s)ds, 
Jr Jm. 



where Z is some continuous stochastic process with K[Z] := zq ^ 0, and 0*(s, ■) is 
the Wick representation of 0(s, ■) = H{2H — l)\s — -\ . Now, we approximate 
ll^e^ + Q\\ — z + Q"^ /2z because Q ~ 0{a), and then evaluate the fractional white 
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noise at 2; + a'^Z{uj): 
z + a^Z) 



= H{2H -l)\z + o?Z - sp^'''-'^ = E{:iE - 1) {{z - + o^zX'^^'''''^ 

we have just took the positive part of the absolute value: it is enough for us examine 
this situation. \i^\o?Z\ = o?zq then 



^ n! [(^ - s) + a2^o] 

and all the terms in the series are of order higher or equal to 2 in a. We just need 
to compare the first term against the deterministic coefficient in the white noise series 
expansion: 

t + o?zq) - t) ~ zo(2// - 2)(i - sf^^-^^o?. (4.32) 

This happens 'coordinate' to 'coordinate' in the fractional white noise decomposition, 
thus we have found 

!^!&)-M^!l±M_.0(„^,. (4.33) 

z Ike^ + Qll 

The first-order equation (4.31) is unaffected by this replacement since these processes 
differ in o? . Finally, we arrived to the desired linear equation: 

W\r,'z)oQ{z) 
J^.Qi-)=9 , (4.34) 

we have set g — a/2 IqUq {g ~ 10~^). 

4.3.1 The stochastic Volterra equation and its solution 

The integral form of equation (4.34) is, 

Q{z) ^Qoz + g / o Q{s) dsds' . (4.35) 

Jo Jo "5 
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Let us set the following initial conditions Q{0) — and Q{0) e — the initial velocity 
is also uncertain. It can be simplified a bit more since 

Jo Jo s 

= / / XloA'')X[o,s')(s) ^"^^'''^ oQ(s) dsds' 
- f ( f Xlo,.)(s')xio,s')(') ds') ^""^^o"'") oQ(s)ds 

= [{z-s)xio,.) ^"^^'''K q{s) ds 
Jr s 

= r ^^^W''{lo^s)oQ{s)ds. 
Jo ^ 

Thus we have a stochastic Volterra equation with (Fredholm) kernel: 



k''{z,s) ■.= g^-^^X[o,.]{s)W^{l^h). 
s 

We will be interested in finding a solution on the (closed) interval < 2; < L. The 
kernel is continous everywhere but s = 0, and 

Wk^'iz, s)\\H,-g < gM X[o,z]{s)s-' \z-s\, (4.36) 

as can be seen from the bound (4.49). 

Now we have to see what are the conditions that make equation (4.35) solvable. It 
should be, if we were able to apply a fixed-point theorem to the above kernel. Therefore, 
proposing as ansatz the usual resolvent for convoluted kernels, that is. 



Kh{z,s)^Y.k'^\z,s), (4.37) 

n=l 
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such that 



Q{z) = Qoz + J Kh{z, s) o (^Qqs ) ds 

= (5oO z+ / KH{z,s)sds 
Jo 



(4.38) 



with the K^^ given inductively by 

K^H^^\z,s)^j KP{z,u)ok^{u,s)du, with n>l, 
K'~i\z,s) = k''{z,s). 



(4.39) 
(4.40) 



It was found by Holden et al. (1996) that this is the unique solution for bounded kernels 
in the distribution Hida space. Their proof is based on the existence of a bound via 
the norm || ■ ||-i -g. The same theorem can also be shown valid in the fractional Hida 
spaces with || ■ Wn-q- But our kernel is unbounded, since the fractional white noise is 
continous and non-zero at s = 0. 

Gripcnbcrg et al. (1990) discuss this type of problematic kernels for normed spaces. 
Defined the space of continous functions / : J — > K with norm ||/||lp(j) = (/j ||/(s) p ds) 
where J C IR is not necessarily compact and X is a Hilbert space — with || • ||. After- 
wards, they introduce a norm for the kernel k: 



LP.p'(J) 




sup / / \\g{z)k{z^s)f{s)\\dzds^ 
liP(j)<i J J J J 

I|9ILp'(j)<1 



p p 



(4.41) 



Then, they proved that a resolvent solution exists whenever this norm is less than one 
(Corollary to Theorem 3.9 in Gripenberg et al., p. 235). This norm has also another 
property, using the Holder inequality the following can be proved^: 



lp.f'(j) < niin< 



Jmz,sw. 



2L 

p 

ds 1 dz 



J \Jj 



\\k{z,s)\\Pdz] ds 



^For a proof see Appendix C. 
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The theorem and property above can be tracked back to the norm in the fractional 
Hida space. Hence, the same hypothesis apphes for this stochastic Predholm kernel 
defined J= (0,L]: 

LP'P'{j),-q < 1 for some q > 0; moreover, 



\\f^"\\Lp,p'iJ),-q < nii 



mm 



J \JJ 



l^(^ljk^{z,sw;,^_^dzy ds 



, (4.42) 



where 11 F| 



LP{J)-q 



{Jj \\F{s)\\^^ _^dsY^^ . Then applying equation (4.36) to the 



bounding condition (4.42) we find*^ 



\\k''\\Lpy(j),-q<9Mp- 



l/p( 



sm np 



i/p' 



<1, 



(4.43) 



since M is a small constant and g <^1. This guarantees the convergence of the proposed 
ansatz. 

Unfortunately, the solution represented as a series of convoluted kernels, eqs. (4.38)— 
(4.40), is useless for calculations. Next, we will prove that a fractional chaos expansion 
exists for the solution. Let us take the second term in the Wick product of equa- 
tion (4.38), it can be written 



X{z)^z + 





oo 



n=l 



s ds 



n=l 



K^^\z, s)s ds 



(4.44) 



because it converges absolutely. The general term in this series can be written, using 
"^Also in Appendix C. 
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definition (4.39), 

r K^^\z,s)sds 
Jo 



l-H 



Tr(n-1) I ^ (52 " Si){Si - 0) 



Sl 



1 

"Z PZ PZ 
Jo J S\J S2 

fZ pz 

]1-H\n 



{S3 - S2){S2 - Si){si - 0) 



S2S1 



J Sn-i Jsn SnSn-1 ■ • ■ Sl 



- (gll-'^r I I I K^" iz, 53) -^v-^ -^v- dss dB^^dBZ 

{Z- Sn){Sn- Sn-l)---{si-0) , r,H 



dBf ■ ■ ■ dBi 



{Si - Si-i) 



~ n 

Jn i=i 

zigz^ll-'^rj f"\s^, ...,si) dB^ . . . dB^^. 
Jwi 



dB^ ■ ■ ■ dB^^ 



(4.45) 



Here we have used the self-similarity property (A. 11) to build the latter adimensional 
integrals, and defined 



/(")(s„,...,si)=(i-son 



1=1 



{Si - Si_i) 



X[Si_l,l)(Si) 



(4.46) 



with So = 0. Now, we build the symmetrized form of the above function, that is, 

/»(s„,...,sO = ^5^/('^H«.„,---,«.J- 



Thus, it induces the following relation 

/ /» (s„, . . . , Sl) dBl ■■■dB^^^j / W (s„, . . . , Sl) dB^ . . . dBl, (4.47) 

JWi_ JWi_ 



because we can rename the each dummy variables of the n\ permutated terms to the 
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normal order. Finally, 



n=l 



(4.48) 



where g = IoQ^z/Iq)^ . This will be nothing else but the fractional chaos expansion 
provided 

oo 

E^'"I/^"^IJ<°° (4-49) 



n=l 



holds. In fact this condition express nothing else that the existence of the variance of 
the process, 



^[X^z)] ^z^ 



1 + 



n=l 



(4.50) 



— we used property (3.84). The search of an upper bound for the succession of 0-norms, 
given that the p'^^ are symmetric, is straightforward: 



1/^ 



n)|2 



/(")(sn, . . . , si)/(")«, . . . , Si)0(s„, 4) • ■ ■ 0(si, si) dSn ' ' ' dsi ds'^- ■ ■ ds[ 
- / / ■■■/ / 4>{sn,s'n)---<l}{si,s[)dsn---dsids'^---ds[, (4.51) 

^0 ^0 Jsn-1 Js'^_^ 



because of definition (4.46) and the fact < Sj — Si-i < Si (idem < — s^_i < s'J 
the last inequality follows. Observing that 

/ / 4>{Sn, O dSnds'^ 

^H{2H-l)j j \sn-sS''~^dsnds'^ 



(1 - Sn-ir + (1 - Ci)'^ - K - s'S" < 1, (4.52) 



we iteratively apply it in (4.51) to find: |/^'*-'|| < 1- Thus, the chaos expansion exists 
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for all 2; < L whenever 

'-©"<^ 

is satisfied. Prom the definition of g and the magnitude of the quantities^ utilized here 
we have: 

L « /^'/=^^(C,2)-V2^^. (4.53) 
So, the condition above is always fulfilled. 



4.4 Ray- light Statistics: a Test Case 

In this section we will use the stochastic ray-equation solution to study the behavior of 
the displacements with respect to the characteristic variables of the system: , /q and 
L. We note that both coordinates of displacement are independent, and they also hold 
the same (non-coupled) differential equation. It is enough to consider a 1-dimensional 
case then. The parameter election (4.17), we have used in our treatment, also defines 
the meaning of the transversal velocities, for they are the angles of deviation. Being 
the velocities continuous we can set. 

Since our solution is dependent of the initial refractive angle 9, its behavior at the 
boundary, e — > 0, should be given. This boundary is just the interface between tur- 
bulent and resting air. Henceforth, we will also model the initial angle as a fractional 
Brownian motion, 

0{e) = c [ X[o,e)i^) dB^ = ci?^(e), (4.54) 

J]St+ 

the constant c is adimensional and measures the strength of the noise. The length 
e works as a kind of correlation distance, as it goes to zero we are examining the 
properties of the interface's short-range correlation. 

Besides, any stochastic process can be put in terms of the spans described in the 
past chapter, and these depend on the construction of stochastic integrals by step 
functions. In any case, even if the former model needs to be corrected — maybe the 
^See pages 3 and 34. 
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interface introduces long-range correlations — the next results are useful; since, they 
are the building blocks for more complex stochastic processes — see p. 101. 

Now, using the chaos expansion (4.48) and the initial conditions given here, the 
solution (4.38) is written : 



Q{z)^0{e)oX{z) 



\ n=l M J 



(4.55) 



Prom the Wick product property (3.71) we see 



E[Q(z)] 

= zcE[B^{e)] E 



oo „ 

1 + / f^''\s^,...,s^)dBZ...dB^^ 

n=l M 



• E[l] = 0. (4.56) 



The evaluation of the variance from experimental data is the most common topic in 
many works related to the optical properties of turbulence because it is directly related 
to the structure constant. Hence, we calculate it setting using property (3.76), 



E[Q\z)] =c^E[(5^(e)oX(z))2] 



E 



iD^.,o,,Mz)r]+E[X\z)] |x[o,)| 



(4.57) 



where we have set X = Y and f = g — X[o,e)- We have already evaluated E[X^(2;)] in 
the latter section. The fractional Malliavin derivative appearing at the right-hand side 
demands elaboration, property (3.74) implies 



D^^ X(z) 



Since the ^-differential is linear we have 



DtX{z)x[o,e){^)ds. 



(4.58) 



DtX{z)^zY,rDl 



f^-\sr.,...,s,)dBl...dBf^ 



(4.59) 



We are going to compute these derivatives now: let us fix n > 2, from the first theorem 
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(3.81) we can commute the stochastic integral and 0-differential, 



f f^^\sr,,...,s,)dBl---dBi 
[ Dt [ f^'^UBl-.-dBl dB^^+ [ f^-UBl...dBlct>{s,s,)ds,. (4.60) 



Now, we recursively commute the operators, the ^-differential and the Wick integral. 
Each time we do so another integral as the last one on the right-hand side of the 
equation above is added. After {n — 1) iterations we reach the innermost integral, thus 
we evaluate 



Dt 



[ /»(s„,...,Si)ciSfJ = / f(-\sn,...,S,)(l>{Sn,s)dSn, 



with the aid of property (3.75). Finally, 



Dt 



H h 



f^'^\s^,...,s,)dB^^---dB^^ 



[ f^^^<l>{s,Sn)ds^dBZ_,---dB^^ + 



f dBl ■ ■ ■ <t>{s, Sk) ds,---dB^^ + ---+ f /» dB^^ . ■ ■ dBfjis, s,)ds. 



— n 



dBZ.r--dB^^, (4.61) 



H 



to arrive to the last equality the symmetry of /*^"^ was employed. Instead, for n = 1 
we just use property (3.75): 



Dl 



F\s^)dB^^ 



f'\si)cl>{s,si)dsi. 



(4.62) 



Afterwards, we can build the fractional Malliavin derivative (3.72) from the series 
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(4.59), 



ds' ds+ 



oo „ 

+ ^(n+l)W 

n=l M 



its second moment is 



/>+^)(s',...,Si) X[o,e)(s)0(s,s') ds'ds 



dBl---dB^\- (4.63) 



E 



2 ~2 

^ 9 



+ 



+ ^ (n + Iff- 



n=l 



f'^^'Hs',-) X[O,e)(s)0(«,s') Cis'rfs 



— ^we used the orthogonal property of these integrals. This series converges, we apply 
the same procedure as before to find a bound for the integrals. What is more, each 
norm appearing in the series is bounded by the zero term. 



f^''^'\s',-) X[0,e){^)<l>{s,s') ds'ds 



< 



[ [{l-s')(p{s,s') ds'ds 
Jo Jo 



= I ^ [1 - e^^+i - (1 - e)'^+^l 



2H + 1 2 J - 4 ^ ^ 



2{2H+1) 

Thus, the existence and uniqueness of (4.63) is guaranteed. Finally, we need the norm 



\xio,e)\l^H{2H-l) / \u-sf''-'duds^e''', 

Jo Jo 
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to calculate the variance of the displacements, 



+ r 



+ 



n=l 



(4.65) 



Now, as the correlation distance goes to zero we recover the initial condition. While 
terms coming from the second moment of X{z) banish (they are all bounded and 
multiplied by e^^), it is not the case with those coming from the fractional derivative. 
We will not go through copious calculations since we are interested in a general outline 
of the solution; thereof, the solution can be expressed as 



E[Q^(.)]=^+ 



n=l 



f'^''+'\s',-)X[0,e)i^mS,s')ds'ds 



(4.66) 



e-»0 



We can estimate a bound for the second term: 



F{~g') = 4j2~g'^^{n + 1] 



n=l 



/^"+'Hs',-)X[O,.)(5)0(5,5') ds'ds 



Mi 



<4f;^^>+l)^ (\) 

n=l ^ ^ 

^ oo 



n=2 



Now, because x/ (1 — = nx^" it is 



n=2 



~9' 



(1 - ~9'y 



(1 - 



and then F(^^) < g"^/ (1 — ^^)^, whenever g < 1. Finally, replacing the values for g, we 
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have 



A 2H+2i2/2.-2H 



(4.67) 



Furthermore, for the range of validity given in the past sections, the contribution of 
the function F is less than 10^^. Thus, the first contribution to Malliavin derivative of 
X completely characterize the variance once the interface's properties are defined. So, 
determining the behavior of the interface is crucial for the present model. 



Conclusions 



We started Chapter 1 making a revision of the up-to-date Passive Scalar Fields prop- 
erties. Also, we have shown the refractive index is among them: this is well-known 
in Atmospheric Optics. Nevertheless, the progress made in Fluid Dynamics on scalar 
turbulence has hardly impacted turbulent propagation. Later on, we compared the 
properties a fair model should comply against those followed by actual optical models. 
Afterwards, we formulated the properties that make the family of isotropic fractional 
Brownian motion a good candidate to simulate the turbulent refractive index: 

• The Structure Function asociatcd to the index n, a scalar field, obeys the power 
law ~ with < if < 1. The value of the (Hurst) parameter H depends 
on the state of the turbulence: ii > 1/2 for highly anisotropic scalar turbulence, 
and H < 1/2, almost always near 1/3, whenever the forces that generate the 
turbulence are not relevant. 

• The Structure Function dependence in r induces a variance corresponding to a 
non-differentiable process. 

• It is assumed a Gaussian process. This is an ad-hoc supposition widely used 
among the literature: it is specially applied when the process plays the role of a 
source in a fluid equation. This approximation is good whenever we are interested 
in the low moments associated to the stochastic process. 

We have proved our proposed model (1.105) fulfills all these conditions. Moreover, 
we obtained its fractal dimension, equation (1.93), matches the estimated by (Con- 
stantin et al., 1991) for passive scalar: dime^^ = 3 — H. Therefore, the exponent H 
determines the state of the turbulence. 
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Finally, we must stress this model give us a local structure function for the refractive 
index — as suggested by Perez and Garavaglia (2001) and some preliminar experimental 
measures. 

On Chapter 2 we have shown under what conditions the wave-equation bring to us 
the paraxial approximation. Then, following Charnotski et al. (1993) we have written 
its Green function using a path integral velocity representation. 

All over this chapter the Markovian approximation is used. It has dominated the 
Atmospheric Optics scenario among the classic models. As it was noted, this model 
discriminates the direction of propagation, from the remaining coordinates. Implying 
a Brownian motion governs the behavior in that direction, that is, 

e(p, z) oc W{z). 

Thus, we used this model to calculate the effects of the turbulence over a system of 
grids (Perez and Garavaglia, 1999). First, we have analized the image formation with 
and without turbulence. We observe the grids arrangement naturaly selects certains 
positions where the visibility is different from zero; that is, the formation of auto- 
images. After the introduction of the turbulence this property remains unchanged. 

On the other hand, the quality of the image is degradated. It depends on the 
geometry of the grids, represented by d and L, as it is shown in figures (2.5) and (2.6). 
In the particular case o? — > cx), the visibility behaves as if the turbulence were absent 
in coincidence with Zavorotny (1988). 

Since the turbulent medium produces a cut-off in Fourier series for the irradiance 
pattern introduces a method to evaluate the structure constant C| as we showed. 

Finally, with the tools exposed in the third chapter we can advance to Chapter 4 
and solve the ray-equation coming from the Geometric Optics in the turbulent case. 

At the introduction to this chapter we have shown substantial differences between 
our model and the Markovian approximation. We also proved that in the markovian 
case it is admisible to commutate derivatives and averages — this is assumed true in 
Optics not caring about the kind of process at hand. Also, this approximation has 
fractal dimension equal to 2|, and thus it is not capable of determine the state of the 
turbulence. There are other models like thise. For example a set of fractal screens 
equispaced has dimension less than 2, and therefore it completely falls out of the 
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foretold range for scalar fields. 

Next, we gave an alternative demonstration to Synge (1937) 's to find ray-equations 
for the (singular) optical lagrangian. The equations for ray light trayectories coming 
from this lagrangian are nonlinear, and then we proceed to linearize them. 

We specifically studied the H > 1/2 problem. The motivations for such a choice are 
various. Prom the mathematical point of view, we were able to define a composition 
of stochastic processes. Afterwards, we have shown the first order ray-equation cor- 
responds to a Stochastic Volterra Equation. Moreover, we have shown that a unique 
analitical solution exists. This solution was expressed as kernel convolutions can be 
rewritten by means of a chaos expansion; thus, turning it into a manageable expression. 

This analysis covers a priori only those cases where average temperature gradients 
are relevant, that is, introduce strong anisotrophies. This behavior is likely to be found 
at the laboratory. Usually, these experiences disregard the process of turbulence mak- 
ing. It is considered that aligning a row of heaters along the ray trajectory (eventually 
using fans) and taking measures at a couple of meters high above them (e.g. Consortini 
et al., 1996) is enough to produce a completly developed turbulence. This asumption 
is at least ingenuous. As we have seen the conditions for isotropy and homogeneity are 
difficult to obtain. First, it must be known for certain the non-existence of a convective 
turbulence; that is, we must observe small Rayleigh numbers for the system (p. 33). 
Also, an inertial tubulence does not necesarly produces an isotropic and homogeneous 
scalar turbulence. As was shown by Villermaux and Gagne (1994), true isotropic and 
homogeneous scalars fields are obtained making the turbulent flow circulate through 
some particular grid arrangement. 

The validity range for our solution contains all the possible distances at the laboratory — 
L <S lO^m. Therefore, our problem is completly determined by the initial conditions; 
our election of the incoming angle as a fBm (4.54) is the right choice given the behavior 
of the scalar quantities. Since this condition is related to short-range correlations we 
should find the constant c depends on the inner scale and the structure constant. 

Afterwards, when we use the solution (4.67) to estimate the variance of a laser 
beam going through the turbulence over a distance L. We obtain: 
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where the correction to this result is of order 0{F) ~ 10~^. Moreover, this term comes 
from the Malliavin derivative. That is, the constant term, of order zero, does not 
contribute to the variance — see (4.65). We must stress that the anisotropy introduced 
by the mean flux should be observed in different constants c at each axis. 

Now, making if + — > 1/2 the displacements variance approachs to Al^^^^L^. This 
is the behavior found by Consortini et al.. It does not correspond to the Kolmogorov 
isotropic model, in accordance to the properties identified at the beginning, but to a 
brownian motion {H = 1/2). That is, given a gaussian process with structure function 
hke in (1.99) the result from Consortini does not hold. This can be achieved when 
Cn \ 1/3; there exists anisotropy o convective turbulence. 

Nevertheless, this result is coherent with the markovian model since the stochastic 
integrals exactly introduce such a dependence with the distance^. The very same 
happens in our case. On the other hand, supposing the extension for H e [1/3, 1/2) 
suggest a similar dependence. It should be changed Iq by L, and thus in that case 
Var Q ~ L^-^^ independently from H. Here we lack the knowledge to estabhsh a value 
for the remaining quantities since the conditions on c are undefined. Although, it is 
clear the power-law difference between this result and the markovian case is relatively 
small. 

These results have been presented in XIII Meeting ON Nonequilibrium Sta- 
tistical Mechanics and Nonlinear Physics (MEDYFINOL'02), December 
9-13, 2002. Another version (Perez, 2002) has been sent to be published. 



Observe 




according to the definition given at the beginning of this section. 
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Fractional Brownian motions 

Before introduce these processes, let us review some basic notions. To build a stochastic 
process a probability space [Q, JF, P] must be provided, where JF) is a measurable 
space with measure P such that = 1 — it is the probability measure. The space 

fl is an abstract space, whose characteristics are irrelevant for the present discussion. 
Now, let {Y, y) be another mensurable space and T a parameter set (e.g., N, R, etc.); 
thus, any given map X : T x fl ^ Y is a. stochastic process if Vt e T 

X-\B) = {u : X{t,uj) e 5} e JP", for any Bey. (A.l) 

Here it will be only necessary to consider y = R and y — B(R) the Borel cr-algebra. 

There is an alternative definition: the canonical representation. We assign to each 
element a; e Q a function X{u;) e R^, where it is defined R^ = {/ : f{t) :T^R}. 
It is called realization of the process. Also, we must provide a cr-algebra so within this 
space the property (A.l) is preserved: the Kolmogorov a -algebra Bi^). It is generated 
by the cylinder sets 

Zt{B) := {/ e R^ : fit) e B}, for B e B(R). 

Finally, from the original probability we can derive the distribution law of X over 
(R^,i3(R^)), 

Fx = P{u; : X{u;) e A}, A e B{R^). 

Therefore, the triad [R^, ;B(R^), Px] constitutes the canonical probability space. The 
original abstract probabihty space [Q, ^, P] is irrelevant if the distribution law of X is 
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given. That is, let us take 

2t„...,u{B, X • • • X S„) = {/ e : ifih), f{tn)) e Si X • • • X e B(R^), 

and thus define the n-dimensional distribution of the process X as Pt^_..._t^(Si x • • • x 
Bn) :— Px(-2^ti,...,t„(-Bi X • • • X Bn)). Conversely, given these finite distributions for all 
n the probability law can be recovered — Kolmogorov's Theorem (Shiryayev, 1984, 
p. 244). 

Henceforth, a Gaussian process can be build from the finite dimensional distri- 
butions, which are normal distributions; that is, \/ti,...,tn £ T the random vector 
{Xti, . . . , Xt^) has distribution 

Pti,...,t„(^l ^ ^ti <Xi + dXi, ...,Xn<Xt„< Xn + dXn) = 



exp --(x-/x)*V-^(x-m) , 



V(27r)'*detV I 2 

where x = (xi, . . . , Xn), /it e and V e x R" is a definite positive matrix. It is 
straightforward to find that {n)i — IE[Xt-] is the mean value at times ti, . . . and 
(V)ij = Cav{Xt^,Xt■) = IE[(Xt. — Hi){Xt. — Hj)] is the associated covariance matrix, 
where E[-] is the average calculated with Fx- Finally, we can formally introduce 
stationarity for processes and its increments. The shift operator is defined (ts o 
f){t) = /(s + t)- A process is called stationary if 

It can be translated in terms of the finite distributions as Pti+s,...,t„+s — Pti,...,t„ for any 
n. Moreover, for Gaussian process this is equivalent to 

Cov(A,„X,J = Cov(A|,,_,.|,Ao), 

for any (ti, . . . On the other hand, a process possess stationary increments if the 
sets r(_|_5 o X — Xi_|_s and o A — A^ has the same distribution. This implies its variance 
has the property 

E [{Xt -Xsf]=E [{X\t-s\ - Xof] . (A.3) 
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= 0, 


almost surely, 


E[B^{s)] 


= 0, 






1 

~ 2 


s + \t\ — \s — t\ 



We are in conditions now to define the 1-dimensional fractional Brownian motion. 
It is a Gaussian process with the following properties (Mandelbrot and Ness, 1968): 

(A.4) 

(A.5) 

(A.6) 

for s, t G M and < H < 1. The exponent H is called Hurst parameter, because it was 
Hurst (1965) who found Nile river's cumulated water flows vary proportional to t^ {t is 
the time) with 1/2 < H < 1. In fact, the family of fBm processes should be separated 
in three subfamilies. When H = 1/2 we recover the standard Brownian motion with 
covariance 

E[B^/\s)B^/\t)] =mm{s,t}:=sAt. (A.7) 
Now, given two dependent Gaussian random variables we have the property 

E{A\B) _ E{AB) 



B E(52) 



(A.8) 



Prom this and the former equation whenever s >t, it is E[B^/'^{s)\B^/^{t)] — B^/^{t). 
This is a martingale, which has no long-memory and its intervals are not correlated. 

On the other hand, the case l/2<iJ<lis the representative case of a long- 
memory process. That is, using equation (A.8) again the conditioned average yields 



1 /Q\2H /q \2H 

E[B''(snB''(t)]^l (£) +l-(^-l) 



B«{t), (A.9) 



i.e., it is not a martingale^. As s grows the conditioned mean behaves 



/ q\ 2H-1 

E[B^{s)\B^{t)]c^H[-) B^{t), 



and diverges at infinity. The long-range dependence is also represented by the diver- 
gence of the series ^^^^ E[B^ {1){B^ (n) - B^{1))] = oo. 

Finally, the case < < 1/2 is left. In the very same way equation (A.9) is 
^If a process X is a martingale it has the property E[X(s)|X(i)] = X{t) for s >t. 
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valid for this range. We observe that E,[B^ {s)\B^ (t)] ~ ^B^{t): on the long-range 
it behaves like a martingale, since it posses short-memory. That is, the correlation of 
the increments is finite < Xli^i ^[B^ {1){B^ (n) — B^{\))\ < oo as the time goes to 
infinity. It is only zero for the Brownian motion, its increments are uncorrelated: it 
has no memory at all. 

The fBm processes have stationary increments. We can evaluate the covariance for 
them. 



1 



I I 1 2iT I I J- J- 1 I ( 4- \ I J. J. 1 2// 

R ~ t^il "T 1^3 ~ ^2| ~ Im ~ ^-21 ~ 1^3 ~ ^'ll 



(A. 10) 



t\\^ ^ and so 



When < ti = ts < ^2 = ^4 we just have E[(5^(t2) - B^{tx)f\ = \t2 
the stationarity is accomplished for the increments. Moreover, if we pick in particular 
— t + h,t3 = t2 = t, and ti = 0, the covariance of the increments according to (A. 10) 

is 

' >0, 



{t + h) 



2H 



t 



2H 



1 



= 0, 

<o, 



\iH> - 
2 

2 

liH <- 
2 



We observe that in the case if > 1/2 consecutive increments tend to have the same 
sign, they are persistent. For the Brownian motion these are as likely to have the same 
sign as the opposite. While in the last case H < 1/2 the increments are more likely to 
have opposite signs, and so we call them anti-persitent. 

As we did with the translation, we define the operator 0q, such that (/ o (f)a)(t) = 
f{oit). We say then the process X is scalar invariant if both X ocj)^ and a^X have the 
same probability distribution for any a and H. For < if < 1 the fractional Brownian 
motion is scalar- invariant: 



B"{as) = a^B^(s), for any a. 



(All) 



where = means they share the same probability law. Usually scalar-invariant processes 
are called self-similar if they have stationary increments. 
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It is worth mentioning that given the change of variable ^ : T ^ T', with $ an 
invective transformation, the redefined stochastic process Zt — X<s>(^t) is also a Gaussian 
process over T with mean = and covariance 

v{s,t)^v'ms),m), (A.12) 

where fi and $ are those defined for the original process X. 

The fractional Brownian motion processes are not differentiable with probability 
1^. This result can be proven from the following lemma found in Cramer and Lead- 
better (1967, §9.4): If a Gaussian process X is differentiable in t with prob. 1 then 
3{d'^v / ds dt){t,t) , v{s,t) is the covariance function of the process. Therefore, the co- 
variance of the fractional Brownian motion (A. 6) implies 

dv ri/[s2^-l-(s-t)2^-l], S>t 

ds"-' ' \H[s^^-^+{t-sf''-\ s<t ^ ' 

Then its diagonal is not derivable and so the second derivative does not exist. 

We have introduced the 1-dimensional fractional Brownian motion process, its n- 
dimensional counterpart can be alternatively constructed through the covariance: 



i=l 

for x,x' e M". 



Fractal Dimension 

The fractal dimension or Hausdorff dimension is defined through the Hausdorff measure 
as follows. 

Given a set F first define a 5-cover as the countable collection of sets {C/j} covering 

^It is said that a process is differentiable with probabihty 1 if P{^n ^ ^} = as n goes to infinity, 
and that a process converges in mean square if E[(^„ — ^)^] as n ^ oo. Both properties are equivalent 
when the process is Gaussian. 
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F, each one with diameter not greater than 5. Henceforth, 

J^s'i^) = |^diam(t/i)^ : {Ui} is a 5-cover of f| , 

where diam([/) = sup^ \x — y\. Then the Hausdorff measure is defined as = 
hm^^o '^'^(-^)- Since this measure is either zero or infinity, the Hausdorff dimension 
of F is univocaUy defined as 

diuiH F = mi{s : J^f\F) = 0} = sup{s : je'{F) = oo}, (A.15) 

and thus, 

{oo if s < dimij F 
(A.16) 
ifs>dim//F 

The direct calculation of this dimension is almost impossible. It is usually done 
through some auxiliar theorems which provides us with upper an lower bounds for the 
Hausdorff dimension. 

In particular lower bounds to the Hausdorff dimension of a set F can be found using 
the potential theory. It is known (Falconer, 1990, Theorem 4.13, p. 64) that given a 
mass distribution /j, on F such that 

fi{dx)fi{dy) 
< oo, 

F - y\ 

it is dimn F < s. The latter integral is known as s-potential. Also, there arc two other 
theorems from Falconer's book we would like to mention here without proof: 
Theorem 7.3: For any sets E cW and F C M™ 




dim/^(E X F) - n < dim^^ + dim^F. (A.17) 



Where dim^F is the upper box-counting dimension defined as, 

diuiBr — hm ; — . 

5^0 - log 6 

where Ns{F) is the smallest number of cubes of side S that cover F. 
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Theorem 8.1: If E,F are Borel subsets ofW^ then 



diuiniE n (F + x)) < max{0, dim.H{E x F) - n} 



(A.18) 



for almost all x e R". 

Theorem 8.2: If E,F C R" he Borel subsets, and let G he a group of transforma- 
tions on R". Then 



for a set of motions a ^ G of positive measure in the following cases: 

(a) . G is the group of similarities and E and F are arbitrary sets. 

(b) . G is the group of rigid motions, E is arbitrary and F is a rectificable curve, 

surface, or manifold. 

(c) . G is the group of rigid motions and E and F are arbitrary, with either <\\va.H E > 



dimnl-E n (t(F)) > dimn E + dimn F -n 



(A.19) 



i(n + l) ordimnF > |(n + l). 
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Markovian Model for the Turbulent Refractive Index 



The markovian model we introduced in Chapter 1 determines a preferred direction 
of propagation, let us say the 2;-axis, and thus the behavior across this direction is 
different from those perpendicular to it. That is, its increments the are independent, 
so they do not have memory of their past. This property, as we mentioned earlier, 
describes a martingale or markovian process. 

Here we will show how the function A in equation (1.101) can be built from the 
original structure function. Let us begin with a locally homogeneous process X having 
as structure function the following: 



Moreover, if we assume it is stationary and Gaussian, as discussed on page 143, its 
correlation function has a spectral representation (Shiryayev (1984), p. 387) 



Because both functions are related by the equation Dx{r) — 2Bx{r) — Bx{0) — -Bx(O) 
we turn the former into 



JR3 

when X is the turbulent refractive index the spectrum is the one discussed earlier in 
Section 1.2.2. 



Dx(r) = (|X(r + r')-X(r')|'). 




(B.l) 




(B.2) 
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Now, taking the inverse transform of (B.l) and using (1.101) we find: 

A{p) = 27r / d^KFx(K,0) 

also, it is 



A{p)^ fdzBx{p,z). 



Besides, when the process is isotropic the spectrum only depends on the absolute value 
of the wavenumber and thus 

poo 

^(p) = 47rM KdK Fx{k,Q) Jo{k\\p\\). (B.3) 
Jo 

Comparing equations (B.2) and (B.3) we define the structure function over the {x,y)- 
plane as 

H{p) = i [^(0) - A{p)] . (B.4) 

TT 

In particular, suppose the power spectra has the form: 

Fx{K,0;z) = ^^^sm^C:{z)\\K\r-\ (B.5) 

then we find 

H{p,z) = S^l±^sm^C^{z)\\pr+\ (B.6) 

^ ^ r[(p + 3)/2]' 2 ^ ^ 

Validity Range of the Path Integral Representation 

To study the validity range of the Feymann's path integral representation we will look 
at the energy flux of a point source radiation through a pupil. If the pupil's transfer 
function is 0(R) then we have 



Poc [ d^RO(R)\G{0,L;R,0)f 



The reciprocity principle implies that this is equivalent to the irradiance / evaluated 
at the point (0, L) provided the initial irradiance distribution function is the same as 
the transfer function. 
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The flux P is a stochastic variable, so we can evaluate its normalized variance or 
scintillation index, that is, 



(7p 



This quantity is an indicator of the type of approximation needed to solve a given 
propagation problem. Since the turbulent refractive index is a Gaussian process the 
mean and free-propagation fluxes coincide, that is, (P) = Pq oc EA;^/47r^L^ with E = 
Jj^2d'^rA{r) the effective area of the pupil. 

Therefore, we can compare the energy flux of the free propagating wave against the 
flux in the turbulent case. Combining equations (2.30), (2.31), and (2.36) evaluated at 
R = we obtain 



in 



where 



and 



4+1 = 



c^VC^(r) // I5V(C)2^S(C)exp[^fc Acvi(C)-V2(C) 



xexp{-$[L,ri(C),r2(C)]}5(2) 



c^Cvi(C) 



dC^2{Q 



, (B.7) 



$ [L, ri(C), r2(C)] = — jj^ {m^M. A - m^2{z), z] - H[v,{z) + v^iz), z] 

with ri{z) — J^d(vi{(), r2{z) — J^d( V2{C) + {1 — z / L) r, and the function H is defined 
as in the former section. Because a more general situation is studied in Chapter 2, we 
have chosen not to give a detailed description for the calculations that lead to equation 
(B.7). 

We observe the strength of the turbulence is measured by the exponential factor in 
the latter equation. Since its arguments have no dimensions, we can show that Up de- 
pends on two dimensionless parameters: the Fresnel number Q = ka'^/L corresponding 
to the pupil effective aperture size a, and q = kpl/L obtained from the spherical wave 
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weak regime 




Figure 2: Tlie graphic displays the weak and strong regime regions. The latter is also 
divided into three subregions: within the region A (f2 -C g <C 1) the scintillation index 
is asymptotically equal to one, in i? (g -C -C q'^^^) is ap = 0{q/il), while in the last 
region (g^/s < < g-i), C, is aj, = 



coherent radius condition D{pq,L) = 1 (Charnotski, 1991, p. 228), where: 

— we used (B.6) for p = 2/3. 

Thus, we define the weak scintillation regime as the set of points {Q, q) where ap 
is asymptotically close to the first term of the Taylor expansion of exp(— $). It is 
found (Dashen, 1979; Zavorotny et al., 1977; Tatarski and Zavorotny, 1980) that this 
condition is reached when q ^ 1 and <^ 1, or g ^ and Q ^ 1. Otherwise, 
the complement to this region corresponds to the strong scintillation regime (Figure 
2). Defined as the region where ap is asymptotically close to the coherent channel 
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expansion. That is produced from two contributions: the main channel expansion 

exp(-$[L,ri(C),r2(C)]) 

= exp 2 dz H[z, v,{z)] I {1 - Q[L, v,{z)M^)] + ■■■}, (B.8) 

where, 



Q[L, ri(^),r2(^)] 

and the additional coherence channel expansion obtained from ri and r2 interchanging 
positions in (B.8). The main idea behind is that the function $ is less than unity in 
one of two regions |ri| ~ po and |r2| ~ po- The coherent channel expansion is thus the 
sum of these two contributions into the scintillation definition 



(7p = M2 + M3 H h A^i + A^2 + ■ ■ ■ 



(B.IO) 



where Mj corresponds to the contribution of the main channel and Ni to that of the 
additional one. By completeness we give here the first terms of these expansions: 



Ni^l f d\ Ca{v) exp [-^(r, z)] 

dz / / d^K F,{k, 0; z) C^(r) 



X sm 



X exp 



1 - 



r • K 



IT 



L 



z z\ L 



dz' H 



(B.ll) 



(B.12) 
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and 



M2 = 27rA;^ J dz / d^K FIk, 0; z) 



A 



X exp 



sm 



(B.13) 



where A is the Fourier transform of the intensity distribution A, is the structure 
function for the index fluctuations, and p{x,y) = mm{x,y}{l — max{x,y}). Further 
terms can be obtained repeating the procedure outhned above. 

This very same procedure can be extended to more complex propagation problems. 
We have seen the scintillation does not only depend on the propagation path but 
also on the aperture size; moreover, the condition Q -C g -C 1, the aperture size 
being larger than the coherence radius, is likely to occur in many situations, even for 
short propagation path, so the only viable tool is the strong-scintillation approach 
(Charnotski, 1996). Thus, the Feymann's path integral approach let us calculate the 
effects of the inhomogeneous media over an irradiance pattern generated by complex 
objects in every possible situation. 



Appendix C 



This Appendix is meant to cover the inequahties shown on page 128. Let || ■ \\H.-q be 
the norm defined in (3.60), F, G, e and J = (0, L]. Using the Holder inequality 
we have 



/ / \\G{s)k''{s,t)F{t)\\H,-,dsdt 

\\G{s)\\h,-J [ \\k''{s,t)ru-,dsY dt 



\H-q,LP{J) 



< 



J \JJ 



Wk^'MWi-qdsY dt 



\F\\H,-q,LP{J)\\G\\H-q,LP'{J)- ■ 1 ) 



for 1/p + 1/p' — l.We apply this very same procedure but beggining with G, thus the 
inequality only has p and p' interchanged. So taking the supremum at both sides yields 



LPP'{J)-q 



< 



J \JJ 



\\k''{t,s)rH^_/s] dt 



J \JJ 



k''it,sW„_dt]- ds 



(C.2) 



Now, we can estimate a bound for the kernel given the above property. Since, we 
know 

||A;^(^, s)\\H,-q < gM X[o,.](s)s-^ |z - s| . 
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Thus, we evaluate: 



t-s\Y 

X[o,t]is]' 1 ds 



dt 



tVdt 



\sm7rp' 



^ (C.3) 
sin np J p 



and 



X[o,t](s)- 



dt 



ds 



(p + l)sP 



ds 



Tip 



V{2p') 



siuTrp' Tip' + i)r(p') y (p + lyip 



(C.4) 



We just compare both terms to reahze that equation (4.43) holds. 
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